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Shervin  Bagheri,  Onofrio  Semeraro,  Espen  Akervik,  Luca  Brandt,  Dan  S.  Henningson 
Linne  Flow  Centre,  KTH  Mechanics,  Stockholm,  Sweden 

1  Introduction 

The  research  efforts  in  the  field  of  flow  control  at  KTH  Mechanics,  led  by  Prof  Dan  Henningson, 
are  partially  supported  through  the  EOARD  Grant  FA  8655-07-1-3053.  Besides  the  project 
leader,  two  graduate  students  are  currently  involved  in  the  project,  Shervin  Bagheri  and  Onofrio 
Semeraro,  and  one  associate  professor,  Luca  Brandt.  Espen  Akervik  who  contributed  to  the 
work  obtained  his  PHD  degree  in  December  2008. 

Control  of  wall-bounded  transitional  and  turbulent  flows  is  the  object  of  several  research 
efforts  owing  to  the  high  potential  benefits.  In  these  fluid-mechanics  systems,  dramatic  effects 
on  global  flow  parameters  may  be  achieved  by  minute  local  perturbations  using  devices  sensing 
and  acting  on  only  small  parts  of  the  flow  with  minute  amount  of  energy.  Such  control  devices 
can  be  used  to  obtain  reduction  of  the  skin-friction  drag,  for  example,  implying  relevant  savings 
of  the  operational  cost  of  commercial  aircrafts  and  cargo  ships.  Our  group  has  been  among 
the  first  to  suggest  and  pursue  the  combination  of  computational  fluid  dynamics  and  control 
theory  [11,  10],  thus  going  past  early  attempts  of  flow  control  based  on  physical  intuition  or  on 
a  trial-and-error  basis. 

In  our  early  work,  a  linear  model-based  feedback  control  approach,  that  minimises  an 
objective  function  which  measures  the  perturbation  energy,  was  formulated  where  the  Orr- 
Sommerfeld  and  Squire  equations  model  the  flow  dynamics.  The  latter  equations  describe  the 
linear  evolution  of  perturbations  evolving  in  a  parallel  base  flow.  The  requirement  implicit 
in  this  formulation  is  the  need  of  complete  state  information.  However,  the  control  problem 
can  be  combined  with  a  state  estimator  to  relax  this  requirement.  The  so  called  Kalman  and 
extended  Kalman  filter  have  been  implemented  in  order  to  reconstruct  the  flow  in  an  optimal 
manner  by  only  considering  continuous  wall  measurements.  Our  studies  have  shown  the  im¬ 
portance  of  physically  relevant  stochastic  models  for  the  estimation  problem  which  turns  out 
to  be  crucial  for  fast  convergence  [9,  6].  Such  stochastic  noise  needs  to  describe  accurately 
enough  the  unmodeled  dynamics,  like  uncertainties  and  nonlinearities.  Based  on  these  models 
the  estimator  is  shown  to  work  for  both  infinitesimal  as  well  as  finite  amplitude  perturbations 
in  direct  numerical  simulations  of  transitional  flows  [5,  13]. 

Motivated  by  the  recent  progresses,  three  aspects  of  the  control  problem  have  been  iden¬ 
tified  as  crucial  in  order  to  apply  feedback  control  in  more  complex  flows  and  achieve  further 
improvement.  They  are  being  addressed  also  thanks  to  the  present  grant  and  are  i)  the  need  to 
go  beyond  the  assumption  of  parallel  base  flows;  ii)  model  reduction  to  significantly  decrease 
the  cost  of  the  estimation  and  thus  of  the  computation  of  the  control  gains;  iii)  the  need  to 
naturally  consider  localized  sensors  and  actuators. 
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1.1  Phase  1 


During  the  first  phase  of  the  project,  feedback  control  and  model  reduction  have  been  applied 
to  the  complex  Ginzburg-Landau  equation,  a  model  problem  for  non-normal  spatial  inhomoge¬ 
neous  flows  which  allows  us  to  consider  globally  unstable  flows  such  as  the  cavity  flow  and  the 
flow  around  a  cylinder  as  well  as  convectively  unstable  flows  such  as  jets  and  boundary  layers. 
The  work  started  on  the  Ginzbug-Landau  model  problem  is  completed  and  the  results  reported 
in  the  review  paper  by  Bagheri  et  al.  [2],  In  addition  matrix-free  methods  are  developed  in 
order  to  compute  the  sets  of  modes  to  be  used  to  design  the  proper  control  algorithm  as  well 
as  to  understand  the  flow  behavior  [4].  These  methods  enable  us  to  compute  two-  and  three- 
dimensional  global  modes  directly  from  the  numerical  solution  of  the  Navier-Stokes  equations. 
In  particular,  Arnoldi  algorithms  can  be  used  to  extract  the  flow  eigenfunctions,  whereas  the 
snapshot  method  is  used  to  retrieve  the  POD  and  the  approximated  balanced  modes,  in  the 
latter  case  using  the  solution  of  both  the  forward  and  adjoint  problem  [15]. 

The  article  by  Bagheri  et  al.  [2]  presents  a  framework  for  the  input-output  analysis,  model 
reduction  and  control  design  for  fluid  dynamical  systems  using  examples  applied  to  the  linear 
complex  Ginzburg-Landau  equation.  Input-output  analysis  generalises  hydrodynamic  stability 
analysis  by  considering  a  finite-time  horizon  over  which  energy  amplification,  driven  by  a  specific 
input  (disturbances/actuator)  and  measured  at  a  specific  output  (sensor),  is  observed.  In  the 
control  design  the  loop  is  closed  between  the  output  and  the  input  through  a  feedback  gain. 
Model  reduction  approximates  the  system  with  a  low-order  model,  making  modern  control 
design  computationally  tractable  for  systems  of  large  dimensions. 

2  Phase  2:  activity  report 

During  the  second  phase  of  the  project,  feedback  control  and  model  reduction  have  been  ap¬ 
plied  to  the  flat-plate  boundary  layer  flow.  The  work  is  completed  and  a  paper  published  in 
the  Journal  of  Fluid  Mechanics  [1].  In  this  article,  the  dynamics  and  control  of  two-dimensional 
disturbances  in  the  spatially  evolving  boundary  layer  on  a  flat-plate  are  investigated  from  an 
input-output  viewpoint.  A  setup  of  spatially  localised  inputs  (external  disturbances  and  ac¬ 
tuators)  and  outputs  (objective  functions  and  sensors)  is  introduced  for  the  control  design  of 
convectively  unstable  flow  configurations.  From  the  linearised  Navier  Stokes  equations  with 
inputs  and  outputs,  controllable,  observable  and  balanced  modes  are  extracted  using  the  snap¬ 
shot  method.  A  balanced  reduced-order  model  is  constructed  and  shown  to  capture  the  input- 
output  behaviour  of  the  linearised  Navier-Stokes  equations.  This  model  is  finally  used  to 
design  a  ^-feedback  controller  to  suppress  the  growth  of  two-dimensional  perturbations  inside 
the  boundary-layer.  This  paper  represents  the  theoretical  foundation  for  the  following  steps 
and  it  is  therefore  enclosed  to  the  report  as  Appendix. 

Prof.  Henningson  gave  the  invited  Julian  D.  Cole  lecture  at  the  AIAA  Theoretical  Fluid 
Mechanics  Conference,  Seattle  June  23-26,  2008  where  he  presented  the  results  on  the  control 
of  boundary-layer  instabilities.  He  was  also  invited  speaker  at  the  European  Fluid  Mechanics 
Conference  in  Manchester  UK  in  September  2008.  An  invited  paper,  based  on  the  invited 
lecture,  is  published  in  AIAA  Journal  [3].  In  this  second  article,  focus  is  particularly  on  the 
development  of  matrix-free  method  for  the  instability  and  control  of  boundary-layer  flows.  It 
is  shown  how  a  numerical  code  solving  the  direct  and  adjoint  Navier-Stokes  equations  can  be 
efficiently  used  to  compute  unstable  modes,  optimal  perturbations  and  balanced  modes  then 
used  to  design  efficient  feedback  control.  The  achievements  in  this  latest  work  which  are  more 
closely  related  to  the  present  project  are  concerned  with  point  iii)  above.  We  have  shown  that 
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Sensors  Actuators 


Figure  1:  Sketch  of  the  input-output  configuration  considered.  The  disturbance  (Bi)  is  modelled 
as  a  localized  TS-wavepacket  located  upstream  at  (xbi,2/Bd2Bi)  =  (20, 1,0).  A  spanwise  row 
(C2  )  of  three  sensors  equally  spaced  along  z  (A2  =  40)  are  used  for  estimation;  the  center  sensor 
is  located  at  (xc2, yc2, ZC2)  =  (150,1,0).  The  actuator  row  (B2  )  has  a  similar  configuration 
with  the  center  actuator  located  at  (xg2,  j/B2,  zs2)  =  (200,1,0).  A  centralized  controller  is 
designed,  i.e.  all  actuators  are  connected  to  sensors  as  shown  by  the  inset  figure.  Further 
downstream  a  region,  that  is  spanned  by  a  number  of  basis  functions  (Ci),  is  used  to  evaluate 
the  disturbance  dynamics  and  thus  acts  as  an  “objective  function”. 


it  is  possible  to  include  in  the  model  wall  sensors  and  actuators,  in  particular  wall-shear  sensors 
and  suction/blowing  holes.  In  view  of  future  practical  realizations,  we  therefore  include  in  the 
simulations  a  more  realistic  actuation. 

The  next  step  towards  more  realistic  realizations  of  feedback  control  were  presented  in  June 
2009  at  the  IUTAM  Symposium  on  laminar /turbulent  transition,  held  in  Stockholm,  Sweden 
[16].  Using  a  number  of  localized  sensors  and  actuators,  a  feedback  controller  is  designed  in 
order  to  reduce  the  growth  of  three-dimensional  disturbances  in  the  flat-plate  boundary  layer. 
A  reduced-order  model  of  the  input-output  system  (composed  of  the  linearized  Navier-Stokes 
equations  including  inputs  and  outputs)  is  computed  by  projection  onto  a  number  of  balanced 
truncation  modes.  It  is  shown  that  a  model  with  50  degrees  of  freedom  captures  the  input- 
output  behavior  of  the  high-dimensional  (n  ~  10')  system.  The  controller  is  based  on  a  classical 
LQG  scheme  with  a  row  of  three  sensors  in  the  spanwise  direction  connected  to  a  row  of  three 
actuators  further  downstream.  The  controller  minimizes  the  perturbation  energy  in  a  spatial 
region  defined  by  a  number  of  (objective)  functions. 

2.1  Input-output  configuration 

The  three-dimensional  input-output  configuration  considered  is  the  extension  of  the  two- 
dimensional  case  studied  in  [1,  3].  We  focus  on  the  dynamics  and  control  of  small  amplitude 
perturbations  about  a  steady  base  flow.  The  main  numerical  tool  is  a  pseudo-spectral  code  that 
provides  solutions  of  the  linearized  Navier-Stokes  equations  and  its  associated  adjoint  equations. 
For  further  details  of  the  code,  boundary  conditions  etc,  we  refer  to  [7,  1]  .  The  computational 
domain  has  the  dimensions  (Lx,  Ly,  Lz)  =  (500,20, 160)  and  a  resolution  of  384  x  81  x  80  grid 
points;  the  fringe  region  starts  at  x  =  400.  The  Reynolds  number  is  Re  =  U0 q5q/v  =  1000, 
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Figure  2:  Streamwise  component  (positive  is  shown  in  black  and  negative  in  gray)  of  the  basis 
functions  Cp*,  with  k  =  1, ...  ,4  given  in  equation  (5). 

where  <5q  is  the  inflow  displacement  thickness. 

The  plant  (shown  schematically  in  Fig.  1),  written  in  an  input-output  form  reads, 


Au  +  Bite  +  B2W 

(1) 

Ciu  +  lu 

(2) 

C2u  +  ag 

(3) 

where  u  is  the  velocity  held,  A  £  Rnxn  is  the  discretized  and  linearized  Navier-Stokes  equations, 
whereas  the  vector  Bi  £  Rnxl  and  the  matrix  B2  £  Rnx3  provide  the  spatial  distributions  of 
the  incoming  disturbance  upstream  and  the  (three)  actuators  downstream.  The  output  signals 
are  extracted  via  the  matrices  Ci  £  Rkxn  and  C2  €  R3xn  that  define  the  spatial  distributions  of 
the  sensors.  The  scalar  quantities  a  and  l  are  penalties  of  measurements  noise  g(t)  and  control 
signal  u(t). 

The  system  (1)  is  stable  since  all  the  eigenvalues  of  A  are  strictly  to  the  left  of  the  imaginary 
axis  on  the  complex  plane.  However,  the  system  is  characterized  by  sensitive  dynamics  as  it 
acts  as  an  amplifier  of  disturbances.  The  upstream  disturbance  consist  of  the  optimal  localized 
initial  condition  computed  by  Monokrousos  et  al.  [12],  that  provides  the  largest  energy  growth 
over  a  given  time. 

The  k  sensors  Ci,  located  further  downstream  are  used  to  define  the  objective  functional 

|Ciu|2  +  l2\u\2  dt.  (4) 

The  aim  of  the  controller  is  to  determine  the  input  signal  u(t),  based  on  noisy  sensor  measure¬ 
ment  v(t)  such  that  the  above  objective  is  minimized.  Note  that  in  this  input-output  framework, 
the  controller  minimizes  the  disturbance  energy  in  a  subspace  of  the  domain,  spanned  by  the 
basis  {Cpi, . . . ,  Cpfc}.  One  choice  of  basis  (the  so-called  output  projection  [15])  are  the  POD 
modes  obtained  from  the  impulse  response  of  all  the  inputs.  This  basis  is  empirical,  i.e.  it 
accurately  represents  the  data  used  to  generate  it.  Using  output  projection  with  few  leading 
POD  modes,  the  controlled  system  shows  significantly  smaller  output  signals  z(t)  compared  to 
the  open-loop  system.  However,  for  three-dimensional  disturbances,  this  does  not  correspond 


4 


Figure  3:  Streamwise  velocity  (positive  is  shown  in  black  and  negative  in  gray)  of  the  first  (top 
row)  and  10th  (bottom  row)  balanced  mode  (left)  and  their  associated  adjoint  mode  (right). 


to  an  actual  reduction  of  the  total  kinetic  energy  of  the  perturbation.  Instead  of  including  a 
very  large  number  of  POD  modes,  alternatively,  a  set  of  spanwise  Fourier  modes  (see  Fig.  2) 
localized  in  the  streamwise  and  wall-normal  directions  can  be  used  as  to  define  a  basis,  of  the 
form 

Ci,fcU  =  jf(0,  exp  (-(x  -  x0)2/al  -  y 2 /a2y,  0)  cos  u  dxdydz.  (5) 

Four  modes  -  from  k  =  1  to  k  =  4  —  localized  around  x  =  300  are  used  in  the  present 
configuration. 

2.2  Model  reduction 

In  order  to  design  a  feedback  controller,  it  is  sufficient  to  capture  the  input-output  (I/O) 
behavior  of  the  system,  rather  than  the  entire  perturbation  dynamics.  For  a  small  number  of 
inputs  and  outputs,  the  I/O  behavior  of  a  linear  system  can  be  described  by  a  reduced-order 
model  obtained  via  the  balanced  truncation  method  [14].  The  reduced  model  retains  the  states 
affected  easily  by  the  inputs  (controllable  states)  and  the  states  that  contribute  the  most  to 
the  outputs  (observable  states).  Essentially,  the  method  amounts  to  an  oblique  projection  of 
the  system  (1),  onto  a  number  of  so-called  balanced  modes  which  can  be  computed  for  high¬ 
dimensional  plants  using  the  snapshots  method  proposed  in  [15]:  snapshots  are  collected  from 
the  impulse  response  of  each  input  via  a  forward  simulation,  and  of  each  output  via  a  simulation 
of  the  adjoint  system  followed  by  one  singular  value  decomposition  (of  the  size  of  number  of 
adjoint  snapshots  times  forward  snapshots). 

Two  balanced  modes  (first  and  10th)  and  their  associated  adjoint  modes  are  shown  in  Fig.  3. 
The  first  balanced  mode  is  nearly  two-dimensional  and  takes  the  shape  of  a  TS  wave-packet 
with  a  large  amplitude  downstream.  This  spatial  structure  is  triggered  with  the  least  energy 
by  the  input  Bi.  Its  corresponding  adjoint  mode  is  essentially  two-dimensional  with  its  largest 
amplitude  upstream.  This  structure,  on  the  other  hand,  generates  the  largest  response  in  the 
sensors  Ci.  The  higher  balanced  modes  (bottom  row  in  Fig.  3)  look  similar  to  the  first  mode 
and  are  characterized  by  different  spatial  wavelengths. 

A  reduced-order  model  of  order  50  is  found  to  capture  the  behavior  of  the  Navier-Stokes 
system  of  order  107.  An  example  of  the  performance  of  the  reduced-order  model  is  shown  in 
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Figure  4:  Impulse  response  from  Bi  — >  Cpi  (right)  and  Bi  — ►  C1.2  (left)  to  a  3D  TS  wavepacket; 
the  solid  line  represents  the  DNS  (n  =  10')  and  the  dotted-line  the  reduced-model  (m  =  50). 

Fig.  4.  With  an  impulse  in  Bi,  a  (optimal)  disturbance  is  introduced  in  the  boundary  layer 
and  grows  as  it  is  convected  downstream.  The  sensor  outputs  Z\(t)  and  zi (t)  extracted  by 
the  sensors  Ci  i  and  Cp2  respectively,  are  shown  with  black  solid  lines:  these  sensors  register 
the  passage  of  a  wave-packet;  the  signal  eventually  decays  to  zero  as  the  disturbance  leaves 
the  computational  box.  In  the  same  figure,  the  output  signals  computed  using  the  reduced- 
order  model  is  shown  (circles)  ,  where  an  impulse  in  the  input  of  the  reduced-order  model  (B i) 
results  in  the  same  response  (extracted  via  the  reduced  sensors  Cit i  and  C\p)  as  the  full  Navier- 
Stokes  system,  albeit  the  significant  order  reduction.  The  approximate  Hankel  singular  values 
(not  shown  here)  decay  rapidely  and  the  leading  singular  values  come  in  pairs  as  observed  in 
previous  studies  [1]. 

2.3  Controller  design 

The  reduced-order  model  can  be  used  to  design  a  controller  of  low-order  that  will  run  “online” , 
next  to  the  numerical  experiments.  Here,  a  classical  Linear-Quadratic-Guassian  (LQG)  (see  e.g. 
[2]  for  introduction  of  control  theoretical  tools  from  a  fluid  mechanics  viewpoint)  is  designed, 
where  all  three  sensors  used  for  estimation  (C2)  are  connected  to  all  three  actuators  B2.  Such 
a  centralized  controller  minimizes  the  energy  of  the  output  signals  (4)  and  more  importantly 
the  resulting  closed-loop  is  guaranteed  to  be  stable.  A  de-centralized  controller  -  when  the 
control  signal  of  each  actuator  is  based  only  on  the  output  from  the  sensor  located  upstream 
and  at  the  same  spanwise  location  -  was  found  both  by  RGA  analysis  [8]  and  by  numerical 
experiments  to  result  in  an  unstable  closed-loop.  This  is  partly  due  to  the  fact  that  localized 
disturbance  introduced  in  the  boundary-layer  spreads  (or  widens)  in  the  spanwise  direction  as 
it  is  convected  downstream,  resulting  in  a  strong  coupling  in  the  spanwise  direction. 

The  performance  of  the  controller  is  shown  in  Fig.  5.  The  level  of  flucuations  (streamwise 
velocity  component  integrated  in  spanwise  and  wall-normal  directions  and  time)  when  the  flow 
is  forced  upstream  with  temporal  white  noise  is  compared  for  three  linearized  DNS;  the  un¬ 
controlled  cases  is  reported  in  solid  black  line  while  the  dash-dotted  and  dashed  lines  show  the 
disturbance  development  when  the  controller  is  active,  i.e.  when  the  measurements  from  the 
three  sensors  (C2)  upstream  are  fed  into  a  controller  that  provide  the  three  actuators  further 
downstream  (B2)  a  control  signal.  The  dashed  line  represents  a  “cheap”  controller  with  l  =  10, 
whereas  the  dashed-dotted  line  is  “expensive”  controller  with  l  =  100  (for  both  controllers 
a  =  0.1).  Near  the  location  (x  =  200)  of  the  actuators  the  growth  of  perturbations  is  trans¬ 
formed  into  a  decay;  further  downstream  the  perturbations  begin  to  grow  again  due  to  the  flow 
instability.  However,  the  overall  disturbance  amplitude  is  reduced. 

In  summary,  a  reduced-order  model  of  order  50  is  able  to  capture  the  input-output  behavior 
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Figure  5:  R.M.S  of  the  streamwise  velocity  component  of  the  uncontrolled  system  (black  line), 
the  cheap  controller  (dashed)  and  expensive  controller  (dotted-dashed  line). 


between  three-dimenstional  disturbances,  actuators,  sensors  and  “objective  functions”.  Using 
this  model,  efficient  control  strategies  can  be  designed  in  order  to  damp  the  growth  of  small- 
amplitude  perturbations  inside  the  boundary-layer.  A  number  of  improvements  are  currently 
under  investigation.  The  spatial  structure  of  the  sensors  and  actuators  will  be  chosen  in  order  to 
reflect  what  actually  can  be  achieved  in  a  practical  experimental  implementation  (for  instance 
with  plasma  actuators). 

3  Conclusions  and  outlook 

The  aim  of  our  EOARD/AFOSR  project  is  to  develop  the  theory  and  tools  to  apply  feedback 
control  in  highly  non-parallel  flow  configurations.  This  can  be  achieved  by  using  global  modes 
to  describe  the  flow  dynamics.  To  design  a  feedback  controller  for  these  more  complex  systems 
it  is  necessary  to  perform  model  reduction.  The  major  difficulty  which  needs  to  be  overcome 
at  this  stage  is  the  derivation  of  a  correct  model  for  the  flow,  the  choice  of  a  successful  control 
algorithms  appearing  as  less  critical. 

During  this  project,  feedback  control  and  model  reduction  have  been  applied  to  the  complex 
Ginzburg-Landau  equation  and  to  the  flat-plate  boundary  layer  flow  [2,  1,  3].  The  results 
indicate  that  in  the  case  of  globally  unstable  flows,  one  or  two  unstable  global  modes  are 
sufficient  to  stabilize  the  closed-loop  system.  However,  for  globally  stable  but  highly  non-normal 
problems  more  modes  are  needed.  In  this  case,  balanced  modes  performed  best  because  of  their 
more  extended  spatial  support:  these  modes  account  for  the  location  of  sensors  and  actuators 
and  do  not  represent  only  the  most  energetic  or  amplified  flow  structures  like  POD  modes  and 
global  modes.  Simulations  of  the  feedback  control  of  perturbations  in  boundary  layers  show 
that  the  control  based  on  a  reduced  order  model  of  the  flow  is  able  to  quench  boundary-layer 
instabilities. 

In  the  near  future  we  aim  at  completing  the  analysis  of  feedback  control  of  three-dimensional 
perturbation  in  boundary  layer  flows  using  localized  sensors  and  actuators.  Recent  developments 
shown  that  it  is  possible  to  include  in  the  model  wall  sensors  and  actuators,  in  particular  wall- 
shear  sensors  and  suction/blowing  holes  [3]  and  that  multiple  arrays  of  actuators  need  to  be 
considered  in  the  design  [16].  The  analysis  presented  here  will  be  extended  to  flows  in  complex 
geometries. 
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The  dynamics  and  control  of  two-dimensional  disturbances  in  the  spatially  evolving 
boundary  layer  on  a  flat  plate  are  investigated  from  an  input-output  viewpoint.  A 
set-up  of  spatially  localized  inputs  (external  disturbances  and  actuators)  and  outputs 
(objective  functions  and  sensors)  is  introduced  for  the  control  design  of  convectively 
unstable  flow  configurations.  From  the  linearized  Navier-Stokes  equations  with 
the  inputs  and  outputs,  controllable,  observable  and  balanced  modes  are  extracted 
using  the  snapshot  method.  A  balanced  reduced-order  model  (ROM)  is  constructed 
and  shown  to  capture  the  input-output  behaviour  of  the  linearized  Navier-Stokes 
equations.  This  model  is  finally  used  to  design  a  Jf2-feedback  controller  to  suppress 
the  growth  of  two-dimensional  perturbations  inside  the  boundary  layer. 


1.  Introduction 

Many  powerful  linear  systems  and  control  theoretical  tools  have  been  out  of  reach 
for  the  fluid  mechanics  community  due  to  the  complexity  of  the  Navier-Stokes 
equations.  Two  elements  that  have  enabled  a  systematic  approach  to  flow  control  are 
the  availability  of  increasingly  powerful  computer  resources  and  the  recent  advances  of 
matrix-free  methods.  In  this  paper,  the  linearized  Navier-Stokes  equations  including 
inputs  and  outputs  are  analysed  using  systematic  tools  from  linear  systems  and 
control  theory.  The  techniques  do  not  rely  on  physical  insight  into  the  specific  flow 
configuration  and  can  in  principle  be  applied  to  any  geometry. 

We  will  focus  on  the  flat-plate  geometry,  which  still  poses  a  computational  challenge. 
The  two-dimensional  Blasius  boundary  layer  is  non-parallel,  i.e.  spatially  evolving, 
and  therefore  has  two  inhomogeneous  spatial  directions.  Many  tools  in  both  stability 
analysis  and  control  theory  rely  on  the  linearized  stability  operator,  which  even 
for  two-dimensional  flows  becomes  very  large  when  discretized.  As  an  example,  a 
moderate  grid  resolution  with  200  points  in  two  directions  leads  to  a  system  matrix 
with  a  memory  demand  of  10  gigabytes,  whereas  storing  a  flow  field  requires  only 
3  megabytes.  It  is  therefore  essential  to  either  approximate  or  develop  algorithms  in 
which  large  matrices  are  avoided,  and  the  storage  demands  are  of  the  order  of  few  flow 
fields.  Matrix-free  methods  employ  the  ‘timestepper  approach'  in  which,  given  a  flow 
field,  a  Navier-Stokes  code  is  used  to  provide  a  field  at  a  later  time.  The  timestepper 
technique  has  become  increasingly  popular  in  stability  analysis,  for  both  computing 
the  largest  transient  growth  (Blackburn,  Barkley  &  Sherwin  2008)  and  performing 
asymptotic  analysis  (Barkley,  Gomes  &  Henderson  2002).  Another  example  of  a 
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Figure  1.  Conceptual  figure  of  the  input-output  configuration  used  for  the  control  of  per¬ 
turbations  in  a  two-dimensional  flat-plate  geometry.  The  computational  domain  Z2  =  (0,  Lx)  x 
(0,  Ly),  shown  by  the  grey  region,  extends  from  x  =  0  to  x  =  1000  with  the  fringe  region  starting 
at  x  =  800.  The  first  input  dSi,  located  at  (xw,  yw)  =  (35,  1),  models  the  initial  receptivity 
phase,  where  disturbances  are  induced  by  free  stream  turbulence,  acoustic  waves  or  wall 
roughness.  The  actuator  provides  a  means  to  manipulate  the  flow,  in  this  case  by  a 
localized  volume  forcing,  and  is  centred  at  (xu,  yu)  =  (400,  1).  Two  sensors  and  ^2  are 
located  at  (x„,  yu)  =  (300,  1)  and  ( xz ,  yz)  =  (750, 1)  respectively.  The  upstream  measurements 
are  used  to  estimate  the  incoming  perturbations,  while  the  downstream  sensor  quantifies  the 
effect  of  the  control.  Note  that  in  this  work  all  the  inputs  and  outputs  are  Gaussian  functions 
given  by  (2.14). 


matrix-free  method  is  the  snapshot  method  introduced  by  Sirovich  (1987),  which 
allows  the  proper  orthogonal  decomposition  (POD)  of  flow  fields  without  solving 
large  eigenvalue  problems. 

The  starting  point  of  modern  optimal  and  robust  control  design,  also  denoted  as 

2-  and  oo-control,  is  an  input-output  formulation  referred  to  as  the  standard  state- 
space  formulation  (Zhou,  Salomon  &  Wu  1999).  The  well-known  stochastic  approach 
to  optimal  control  referred  to  as  linear  quadratic  Gaussian  (LQG)  is  an  example 
of  an  controller.  In  this  work,  we  consider  three  inputs  and  two  outputs;  the 
inputs  represent  external  disturbances,  measurement  noise  and  the  actuator,  whereas 
the  outputs  represent  measurements  for  estimation  and  of  the  objective  functional  to 
be  minimized  (see  figure  1).  The  control  problem  is  to  supply  the  actuator  with  an 
optimal  signal  based  on  the  measurements  taken  from  the  first  sensor,  such  that  the 
effect  of  external  disturbances  and  measurement  noise  on  the  disturbance  energy  is 
minimized  at  the  location  of  the  second  sensor.  Given  the  physical  distribution  of 
the  inputs  and  outputs,  the  control  design  process  amounts  to  the  determination  of 
input  signals  when  output  signals  are  given.  Therefore,  for  successful  control  design 
it  is  sufficient  to  capture  only  a  fraction  of  the  dynamics,  namely  the  relationships 
between  input  and  output  signals. 

The  aim  of  this  study  is  to  build  a  model  of  low  dimension  that  captures  the  input- 
output  behaviour  of  the  flat-plate  boundary  layer  and  use  this  model  for  optimal 
feedback  control  design.  With  the  help  of  the  adjoint  Navier-Stokes  equations  two 
fundamental  dynamical  structures  are  identified:  (i)  the  flow  structures  that  are 
influenced  by  the  inputs;  (ii)  the  flow  structures  that  the  outputs  are  sensitive  to. 
These  controllable  and  observable  structures  determine  the  input-output  behaviour 
completely  for  linear  systems.  It  is  well  known  in  systems  theory  that  these  two 
sets  of  modes  can  be  balanced  and  represented  by  one  set  of  modes,  called  the 
balanced  modes.  In  this  way,  the  flow  structures  that  capture  most  of  the  input- 
output  behaviour  are  extracted  and  used  as  the  projection  basis  for  model  reduction. 
The  method  employed  in  this  work  to  compute  the  balanced  modes  is  called 
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snapshot-based  balanced  truncation  (Willcox  &  Peraire  2002;  Rowley  2005).  This 
method  has  been  applied  to  the  channel  flow  (Ilak  &  Rowley  2008)  and  the  flow 
around  a  pitching  airfoil  (Ahuja  et  al.  2007).  Unlike  previous  work,  we  do  not  combine 
the  snapshot-based  balanced  truncation  with  an  output  projection  approach  in  order 
to  describe  the  flow  dynamics.  Our  control  design  and  performance  evaluation  is  based 
on  input  and  output  signals  rather  than  on  the  space-time  evolution  of  the  entire  flow. 

Previous  work  in  flow  control  involving  model  reduction  and  control  design  has 
typically  relied  on  physical  insight  into  the  specific  flow  situation  rather  than  a 
systematic  approach  detached  from  the  application  (see  Kim  &  Bewley  2007,  for  a 
recent  review).  For  parallel  flows,  for  instance  it  is  possible  to  decouple  the  linear 
equations  in  Fourier  space.  Control,  estimation  and  other  types  of  optimization  can 
then  be  performed  independently  for  each  wavenumber  and  transformed  back  to 
physical  space.  This  approach  has  been  adopted  for  channel  flow  in  Flogberg,  Bewley 
&  Henningson  (2003)  and  even  extended  to  weakly  nonparallel  flows  by  Chevalier 
et  al.  (2007a).  Another  example  is  the  projection  of  the  linearized  Navier-Stokes 
equations  on  a  set  of  modes  such  as  global  eigenmodes  of  the  stability  operator  or 
POD  modes.  Although  these  methods  have  been  applied  with  considerable  success  to 
various  flows  (Gillies  1998;  Akervik  et  al.  2007)  their  success  is  strongly  dependent 
on  the  dynamics  of  the  specific  flow  situation.  For  many  open  shear  flows  the 
global  eigenmodes  and  their  associated  adjoint  modes  can  become  widely  separated 
in  the  stream  wise  direction  (Chomaz  2005)  and  gradually  move  away  from  the 
locations  of  the  inputs  and  outputs  (Lauga  &  Bewley  2003).  As  a  consequence 
controllability  and  observability  of  the  global  eigenmodes  is  gradually  diminished.  If 
controllability/observability  is  lost  for  any  unstable  eigenmode,  no  control  scheme 
will  be  able  to  stabilize  the  system.  The  POD  basis  also  has  limitations  for  describing 
the  input-output  behaviour.  Although  it  is  optimal  for  capturing  the  energy  of 
the  response  to  an  input,  it  does  not  always  capture  the  input  itself  and  takes  no 
consideration  of  the  output.  Flowever,  examples  of  successful  adaptations  of  POD 
modes  can  be  found  in  Noack  et  al.  (2003)  and  Siegel  et  al.  (2008)  for  the  globally 
unstable  flow  past  a  circular  cylinder. 

The  paper  is  organized  as  follows:  We  start  with  describing  the  flow  domain,  the 
inputs,  the  outputs  and  the  control  problem  in  §2.  In  this  section  the  mathematical 
framework  is  presented  with  evolution,  controllability  and  observability  operators 
and  their  associated  adjoint  operators.  These  operators  are  used  to  introduce  the 
Gramians  and  balanced  modes  in  §3,  where  we  also  investigate  the  input-output 
behaviour  of  our  linear  system  and  discuss  the  controllable,  observable  and  balanced 
modes.  In  §  4  the  impulse  and  harmonic  response  of  the  balanced  reduced-order  model 
(ROM)  are  compared  to  the  full  Navier-Stokes  equations,  and  the  model  reduction 
error  is  quantified.  Section  5  deals  with  the  control  design.  We  briefly  introduce  the 
Jf2  framework  and  evaluate  the  closed-loop  performance.  Concluding  remarks  and 
a  summary  of  the  presented  material  are  offered  in  the  last  section.  Finally,  in  the 
appendices  we  derive  the  adjoint  operators  and  describe  the  snapshot  method,  the 
solution  of  the  JC2  problem  and  our  timestepper. 

2.  Problem  formulation 

2.1.  Governing  equations 

We  consider  the  linear  spatio-temporal  evolution  of  two-dimensional  disturbances  in 
a  viscous,  incompressible  flow  over  a  flat  plate.  The  geometry  of  the  problem  and 
the  physical  domain  D  =  (0,  Lx)  x  (0,  Lv)  are  shown  in  figure  1.  The  disturbance 
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behaviour  is  governed  by  the  Navier-Stokes  equations  linearized  about  a  spatially 
evolving  zero-pressure-gradient  boundary  layer: 

M  —  — (£7  •  V)m  —  (u  •  V)t/  —  Vp  +  Re  *V2m  +  k(x)u,  (2.1  a) 

at 

0  =  V  •  u,  (2.1b) 

u  =  u0  at  t  =  0.  (2.1c) 

The  disturbance  velocity  and  pressure  field  at  position  x  =  (x,  y)  and  time  t  are 
represented  by  u(x ,  t)  =  (u,  v)T  and  p(x,t),  respectively.  The  divergence  operator  is 
denoted  by  V  =  (3X,  3v)r.  The  Reynolds  number  is  defined  as  Re  =  U^S^/v,  where  U& 
is  the  free  stream  velocity  and  8 J  the  displacement  thickness  at  the  computational 
inflow  jc0  =  0.  All  the  simulations  were  performed  at  Re  =  1000  which  corresponds  to 
a  distance  of  341<5q  from  the  leading  edge  to  the  inlet  of  the  computational  domain. 
The  base  flow  U  =  (U,  V)T(x,  y)  is  a  solution  to  the  steady,  nonlinear  Navier-Stokes 
equations. 

The  term  X(x)  is  used  to  enforce  periodicity  of  the  physical  flow  in  the  stream  wise 
direction,  so  that  a  spectral  Fourier  expansion  technique  can  be  employed  for  our 
numerical  solution.  This  function  is  non-zero  only  in  a  fringe  region  at  the  end  of 
the  domain  (see  figure  1)  in  which  it  forces  the  outgoing  perturbation  amplitude  to 
zero  (see  Appendix  C  and  Nordstrom  Nordin  &  Henningson  1999  for  further  details). 
As  discussed  in  Akervik  et  al.  (2008)  the  growth  rate  of  individual  eigenvalues  in 
the  spectrum  of  the  linearized  Navier-Stokes  equations  depends  on  the  outflow 
boundary  conditions.  However,  the  perturbation  dynamics  remain  unaltered  for 
different  boundary  conditions,  including  the  fringe  method. 

The  solutions  to  (2.1)  satisfy  no-slip  condition  at  the  plate  and  vanish  at  the  upper 
boundary  Lv  =  308q  which  is  chosen  to  be  well  outside  the  boundary  layer.  The 
boundary  conditions  hence  are 


m(0,  y)  =  u(Lx,  y),  (2.2 a) 

u(x,  0)  =  u(x,  Ly)  =  0.  (2.2 b) 


2.2.  Standard  state-space  formulation  and  the  3^2  problem 
The  Navier-Stokes  equations  may  be  written  in  the  ‘standard  state-space  form’  (Zhou, 
Doyle  &  Glover  2002)  useful  for  applying  tools  from  systems  theory  and  for  /  JfA 

control  design. 

In  the  state-space  framework,  any  divergence-free,  smooth  disturbance  field  u(x) 
that  satisfies  the  boundary  conditions  (2.2)  is  an  element  of  the  (Hilbert)  state  space 

X  =  (m(jc)  G  L2(D)  |  V  •  u(x)  =  0,  u(0,  y)  =  u(Lx,  y),  u(x,  0)  =  m(x,  Ly)  =  0}.  (2.3) 

A  state  is  a  velocity  field  u(x,t)  at  time  t  or  equivalently  a  point  on  a  trajectory  in 
X.  Let  us  introduce  a  bounded  linear  solution  operator  &~(t)  :  X  — >  X  for  the  state 
variable  u  as 


u(x,t  +  s)  =  3T(t)u(x,s).  (2.4) 

Given  a  perturbation  field  at  time  5,  t7~(t)  provides  the  velocity  field  at  a  later  time 
t  +  s  by  solving  (2.1)  with  «0  =  u(x,  s)  and  (2.2).  The  operator  satisfies  the  properties 
3~(t  +  $)  =  3T(t)3F(s),  2T (0)  =  /  and  can  be  considered  a  semi-group  (see,  e.g.  Pazy 
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1983;  Trefethen  &  Embree  2005)  of  the  form  exp(,a?t)f  with  the  infinitesimal  generator 

3T(St)u-u  (25^ 


sslu  =  lim 
St-*  0 


St 


The  linearized  Navier-Stokes  equations  (2.1)  with  boundary  conditions  (2.2)  can  be 
cast  as  an  initial-value  problem  in  state-space  form 


ii  =  s/u, 
u  =  u0 


at  t  =  0. 


(2.6) 

(2.7) 


Note  that  the  action  of  si  on  u  corresponds  to  evaluating  the  right-hand  side  of 
the  Navier-Stokes  equations  and  enforcing  the  boundary  conditions.  The  pressure 
term  can  be  obtained  from  the  velocity  field  by  solving  a  Poisson  equation  (Kreiss, 
Lundbladh  &  Henningson  1993).  Alternatively,  the  state-space  form  can  be  obtained 
by  defining  a  projection  operator  that  projects  the  equations  onto  X  (Chorin  & 
Marsden  1990;  Bewley,  Temam  &  Ziane  2000). 

In  this  work,  the  action  of  the  operator  si  is  approximated  numerically: 
&~(t)u(x,s)  is  obtained  by  solving  the  partial  differential  equation  (2.1)  using  a 
timestepper,  i.e.  a  Navier-Stokes  solver  (Barkley  et  al.  2002),  with  m(jc,s)  as  initial 
condition.  In  its  simplest  form,  a  timestepper  sets  up  a  grid  in  space  and  time  and 
computes  approximate  solutions  on  this  grid  by  marching  in  time.  This  approach  is 
computationally  feasible  also  for  very  large  systems,  since  matrices  are  not  stored. 
The  timestepper  used  and  the  corresponding  numerical  method  are  described  in 
Appendix  C. 

We  introduce  the  forcing  f(x,  t),  which  is  also  referred  to  as  the  input.  The  forcing 
/  is  decomposed  into  external  disturbances  3Sy  w  and  a  control  i.e. 

/  =  fjw  +  f2u,  (2.8) 

where  the  input  signals  w(t),  u(t)  are  functions  of  time  and  are  bounded  linear 

mappings  from  IR  — *•  X.  The  first  mapping,  represents  the  spatial  distribution  of 
the  sources  of  external  disturbances  acting  on  the  flow  (see  figure  1).  In  our  model,  the 
input  forcing  3Sy  is  located  at  the  upstream  end  of  the  domain  to  model  the  upstream 
receptivity  phase,  when  disturbances  are  introduced  into  the  boundary  layer  by,  e.g. 
roughness  and  free  stream  perturbations.  The  actuator  used  for  control  is  defined  by 
the  mapping  0S2,  which  represents  a  localized  volume  force,  mimicking  blowing  and 
suction  at  the  wall.  Finally,  u(t)  represents  the  control  signal  we  wish  to  apply  and  is 
based  on  the  sensor  measurements. 

Information  about  the  disturbance  behaviour  is  given  by  two  outputs : 

Z  =  ^yU  +  lu,  (2.9) 

v  =  u  +  ag,  (2.10) 


where  the  output  signals  z,  v  are  functions  of  time  and  (4\,  (42  are  bounded  linear 
functionals  from  X  — *  IR.  The  sensor  defined  by  'Hy  is  located  far  downstream,  and 
it  is  used  to  evaluate  the  level  of  the  disturbance  amplitude.  Therefore  it  reveals 
whether  the  ‘objective’  of  our  control  has  been  met.  In  particular,  the  objective  is  to 
find  a  control  signal  u(t)  such  that  the  perturbation  energy  in  the  flow  is  minimized 
downstream  at  the  location  defined  by  'Hy.  To  design  an  efficient  controller,  however, 


f  This  is  only  true  if  si  is  a  bounded  operator.  In  general,  for  a  closed  operator  si  with  a  dense 
domain,  the  relation  between  and  si  is  i~(t)u  =  lim,,^^/  —  t/nsifu  (see,  e.g.  Pazy  1983). 
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the  energy  input  expended  in  the  actuation  should  be  limited;  thus,  the  control  effort 
is  penalized  with  a  scalar  1.  For  large  values  of  /  the  control  effort  is  considered  to 
be  expensive,  whereas  small  values  indicate  cheap  control.  This  results  in  an  objective 
functional  of  the  form 


||Z||  =  ||<?m||+/||u||  (2.11) 

and  explains  why  the  control  signal  is  added  to  the  sensor  signal  when  dehning  the 
output  signal  z.  The  norms  in  (2.11)  are  associated  with  the  inner  products  defined  in 
the  next  section.  In  the  definition  of  z  we  have  assumed  (lu,  Wiu)  =0,  so  that  there  is 
no  cross-weighting  between  the  flow  energy  and  control  input  (Zhou  et  al.  1999). 

The  second  output  signal  v{t)  is  the  measurement  signal  extracted  from  the  sensor 
e&2-  This  signal  is  the  only  information  delivered  to  the  controller  in  order  to  provide 
a  control  signal  such  that  the  above  objective  is  met.  The  additional  term  g(t) 
accounts  for  noise  contaminating  the  measurements.  This  term  can  be  considered 
as  a  third  forcing,  but  rather  than  forcing  the  Navier-Stokes  equations  it  forces  the 
measurements.  Large  values  of  the  scalar  a  indicate  high  level  of  noise  corruption  in 
the  output  signal,  whereas  for  low  values  of  a  the  measurement  v  reflects  information 
about  the  flow  field  with  high  fidelity. 

The  choice  of  the  relative  position  of  the  sensor  # 2  and  actuator  used  in 
the  control  design  process  and  reported  in  figure  1  is  based  on  the  knowledge 
of  the  behaviour  of  boundary  layer  instabilities.  For  convectively  unstable  flows, 
disturbances  eventually  leave  the  control  domain;  therefore  there  exists  only  a  window 
of  opportunity  in  time  to  reduce  the  growth  of  these  disturbances,  while  they  are 
convected  downstream.  As  a  consequence,  the  sensor  placement  in  the  streamwise 
direction  is  a  trade-off.  For  a  good  estimation  performance  it  should  be  placed 
downstream,  so  that  the  disturbance  energy  has  amplified,  but  for  a  good  controller 
performance,  it  should  be  placed  upstream  to  provide  the  actuator  an  estimate  of  the 
flow  dynamics  as  soon  as  possible.  Similarly,  there  is  a  trade-off  when  choosing  the 
location  of  the  actuator,  since  its  effects  on  the  disturbance  behaviour  is  limited  to 
the  nearby  region.  It  is  rather  inefficient  to  place  it  either  far  downstream,  where  the 
disturbances  have  already  experienced  a  substantial  growth,  or  far  upstream,  where 
the  disturbances  will  again  have  the  opportunity  to  grow. 

A  completely  different  choice  of  sensor  and  actuator  placement  is  appropriate  in 
the  case  of  globally  unstable  flows  (see  Bagheri  et  al.  2008),  when  the  whole  flow 
beats  at  a  specific  frequency.  Since  the  disturbances  never  leave  the  laboratory  frame, 
one  can  place  the  measurement  sensor  at  the  place  at  which  the  disturbance  energy  is 
the  largest  and  the  actuator  at  the  place  at  which  the  sensitivity  of  the  disturbances 
to  forcing  is  the  largest.  In  many  open  shear  flows  these  locations  are,  respectively, 
far  downstream,  where  the  global  eigenmodes  of  the  linearized  operator  are  located, 
and  far  upstream,  where  the  adjoint  global  eigenmodes  reveal  high  flow  sensitivity. 
See,  for  example  Akervik  et  al.  (2007)  for  the  case  in  which  feedback  control  applied 
to  flow  separation  over  a  long  shallow  cavity. 

The  Navier-Stokes  equations  (2.1)  with  input  vector  f  =  (w.g.u)T  as  an  element 
in  the  input  space  IU  =  R3  and  output  vector  y  =  {z ,  v)T  as  an  element  in  the  output 
space  Y  =  R2  may  now  be  written  in  the  standard  state-space  form 


ii  =  s/u  +  %f, 
y  =  (€u  +  <2)f, 

U  =  uq  at  t  =  0, 


(2.12u) 

(2.126) 

(2.12c) 
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Figure  2.  The  operators  used  to  examine  the  system  input-output  behaviour.  The 
controllability  operator  Y£  c  relates  past  inputs  to  the  present  state,  while  the  observability 
mapping  Y£a  relates  the  present  state  to  the  future  outputs.  Their  combined  action  is  expressed 
by  the  Hankel  operator  . 


where  si  has  been  defined  in  (2.5).  Furthermore,  we  have  <€  =  (^i,  ,  28  = 

J2)  and 

fo  o  A 

@=\  .  (2.13) 

\0  a  OJ 

The  system  (2.12)  is  asymptotically  stable;  i.e.  in  the  global  framework  all  the 
eigenmodes  of  the  linearized  Navier-Stokes  system  for  a  spatial  boundary  layer 
represent  perturbations  decaying  in  time. 

Finally,  we  define  the  spatial  distribution  of  the  sensors  and  actuators  introduced 
above.  Here,  the  input  and  output  operators  are  modelled  with  the  Gaussian  function 
h(x,xo),  defined  as 


where 


h(x\x  0) 


/  °xYy 

y-wn 


exp  ( -y; 1 


(2.14) 


X  —  XO 


Yy  = 


y-yo 

dy 


(2.15) 


The  scalar  quantities  ox  =  4,  crv  =  1  /4  and  x0,  yo  (the  latter  two  being  given  in  the 
caption  of  figure  1)  determine,  respectively,  the  size  and  location  of  the  inputs  and 
outputs.  They  are  all  of  the  same  size  but  located  at  different  streamwise  locations, 
as  shown  schematically  in  figure  1.  With  these  definitions  we  have 


28  =  (/i(x;xw),  0,  li{x\xu)) 


(2.16) 


and 


(h(x,xz)Tu 
h{x,xv)Tu 


dx  dy. 


(2.17) 


The  particular  shape  of  sensor  and  actuators  implies  that  the  inputs  amount  to 
localized  volume  forcing,  whereas  the  flow  measurements  are  obtained  by  averaging 
the  velocity  field  over  small  domains,  using  the  Gaussian  function  as  weights. 


2.3.  Controllability  and  observability  operators 

When  performing  model  reduction  for  control  design,  one  wishes  to  retain  the 
relationship  between  the  inputs  and  the  outputs  in  the  low-order  system.  Following 
linear  systems  theory,  the  properties  of  the  input-output  system  (2.12)  can  be  described 
by  the  two  operators  introduced  in  this  section.  In  the  framework  presented  next  we 
assume  that  sufficient  regularity  exists,  so  that  all  operators  are  bounded  in  the  chosen 
metrics.  See  figure  2  and  table  1  for  an  overview  of  the  operators. 
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Definition 


Adjoint  operator 

t7~(t)u(s)  =  u(t  +  s)  =  n(s  —  t) 

ro 


Operator  Mapping 

Evolution  X  — >  X 

Controllability  U((-oo,  0] )  ->  X  &cf(t)  =  #~(-t)&f(t)  At  £Tcu0  =  &'f(-t)u0 

Observability  X ->•  Y([0,  oo))  £Co(t)u0=^(t)u0  JSf*(f)/=  /"  F’(tyg'iAt 

Hankel  U((-oo,  0])  ->  Y( [0,  oo))  y  =  Seo(t)Secf  f=£Tc(t)£roy 

Table  1.  The  linear  operators  used  in  this  work.  See  Appendix  A  for  further  details  and 
derivations  of  the  adjoint  operators. 


The  operators  needed  to  describe  the  input-output  behaviour  can  be  related  to  the 
formal  solution  of  the  system  of  (2.12),  which  is 

y(t)  =  c&^~(t)u0  +  <€  f  2T(t  —  r)^f(r)  dr  +  3>f{t).  (2.18) 

Jo 

In  this  expression,  we  identify  the  first  term  on  the  right-hand  side  with  the 
homogeneous  solution  and  the  second  term  with  the  particular  solution  stemming 
from  the  forcing  f.  Note  that  in  our  case  the  forcing  term  tffl  is  an  element  in  X; 
i.e.  it  is  divergence-free  and  satisfies  the  boundary  conditions.  For  a  more  general 
forcing  /,  only  the  divergence-free  part  of  the  forcing  /  will  affect  the  output  signal. 
The  difference  /  —  /  can  be  written  as  the  gradient  of  a  scalar  and  thus  will  only 
modify  the  pressure  (Bewley  et  al.  2000).  The  third  part  of  (2.18)  relates  the  input 
to  the  output  through  the  matrix  S>  without  any  operators  involved.  Without  loss  of 
generality,  we  will  neglect  this  term  for  now  and  consider  it  again  in  §  5  for  control 
design. 

In  systems  theory,  the  quantitative  investigation  of  the  input-output  properties  of 
a  linear  system  is  commonly  performed  through  the  mappings  sketched  in  figure  2. 
We  begin  by  introducing  the  controllability  operator  J5?c  :  TU(  ( —  co,  0])  ->  X: 

u0  =  £'cf{t)  =  f  3T(- r)^/(r)dr.  (2.19) 

J  —  OO 

This  operator  describes  the  mapping  of  any  input  f{t)  with  t  e  (— oo,  0]  onto  the 
state  vector  u  at  the  reference  time  t  =  0.  The  input  space  UJ((— oo,  0])  contains  input 
trajectories  in  the  past  time  t  e  (— oo,  0],  The  associated  inner  product  is  given  in  §  A.l 
of  Appendix  A.  The  action  of  J5? c  can  be  numerically  computed  by  a  timestepper. 
It  amounts  to  solving  the  linearized  Navier-Stokes  equations  for  the  velocity  field  u 
with  forcing  term  f(t)  and  zero  initial  conditions. 

The  observability  operator  T?0  :  X  — >  Y([0,  oo))  is  defined  as 

y(t)  =  J5? o{t)u0  =  ^^~(t)u0.  (2.20) 

This  operator  describes  the  mapping  of  any  initial  velocity  field  u0  to  the  output 
signal  y(t)  with  t  ^  0.  The  output  space  Y([0,  oo))  contains  output  trajectories  in  the 
future  time  t  6  [0,  oo).  The  action  of  =S?0(f)  can  also  be  numerically  computed  and 
it  amounts  to  extracting  the  output  signal  while  solving  the  linearized  Navier-Stokes 
equations  with  the  initial  condition  u0  at  the  reference  time  t  =  0  and  zero  forcing. 

A  direct  mapping  between  input  and  output  can  be  obtained  as  the  combination 
of  the  operators  just  introduced  (see  figure  2) : 

y(t)  =  J£„Secf(t)  =  f  V2~(t-  r)#f(r)dr.  (2.21) 

J  —  OO 
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This  expression  can  be  interpreted  as  a  mapping  from  past  inputs  to  future  outputs. 
It  can  be  shown  that  (2.21)  is  the  formal  solution  for  a  system  which  is  forced  by  f(t ) 
in  the  time  interval  t  G  (—00,  0],  resulting  in  the  flow  field  u0  at  t  =  0.  The  output  is 
extracted  for  t  ^  0,  corresponding  to  the  signal  y(t)  produced  by  the  initial  condition 
u0.  The  expression  (2.21)  is  also  the  starting  point  for  the  input-output  analysis, 
leading  to  systematically  finding  reduced-order  approximations.  The  mapping  from 
the  inputs  to  the  outputs  given  by  (2.21)  in  terms  of  and  is  called  Hankel 
operator  :  UJ((— 00,  0])  — >  Y([0,  00)),  i.e. 

m  =  &0(t)&cf(t)  =  (Jff)(0.  (2.22) 

We  have  two  different  representations  of  the  input-output  behaviour  of  the  flow 
system:  (i)  the  state-space  representation  (2.12)  and  (ii)  the  Hankel  operator 
defined  in  (2.22).  Note  that  in  the  latter  case  it  is  assumed  that  the  inputs  and  outputs 
are  not  active  at  the  same  time. 

2.4.  Adjoint  equations  and  operators 

Before  issues  related  to  controllability,  observability  and  model  reduction  can  be 
considered  the  adjoint  linear  operators  corresponding  to  {2T ,  jS?c,  j£?0)  must  be 
introduced.  The  adjoint  variables  provide  information  about  how  variations  in 
the  velocity  field  affect  the  system  output.  We  show  that  the  adjoint  operators 
can  be  associated  with  the  adjoint  linearized  Navier-Stokes  equations  in  state- 
space  form,  where  the  role  of  the  inputs  and  outputs  is  reversed.  The  operators 
IflT* ,  jS?*,  jS?*,  <$* ,  38*)  and  adjoint  Navier-Stokes  equations  are  derived  in  Appendix  A. 
The  inner  products  associated  with  the  Hilbert  spaces  X,  UJ,  X,  HJ(( — 00,  0])  and 
Y([0,  00))  are  also  given  in  §A.l  of  Appendix  A.  The  adjoint  of  the  linearized 
Navier-Stokes  equations  (2.1)  associated  with  inner  product  (A 2a)  is 

— p  =  (U  •  V)/>  —  iflU)T p  +  Vcr  +  Re~lV2  +  k(x)p{ jc),  (2.23 a) 

0  =  V‘P,  (2.23  b) 

p  =  pT  at  t  =  T,  (2.23c) 

This  system  of  equations  describes  the  evolution  of  adjoint  flow  field  p(x,  t)  =  (u* ,  v*)T 
backwards  in  time.  The  term  a  denotes  the  adjoint  pressure  field.  The  boundary 
conditions  for  p  are  given  in  §  A.  1  of  Appendix  A. 

The  evolution  operator  associated  with  (2.23)  is 

p(x,  s  —  t)  =  3~*(t)p{x,  5),  (2.24) 

so  that  given  an  adjoint  field  at  time  s  the  adjoint  evolution  operator  provides 
a  solution  at  an  earlier  time  s  —  t.  Again,  the  above  operator  is  approximated 
numerically  using  a  timestepper  solving  (2.23).  In  §A.l  of  Appendix  A  it  is  shown 
that  3T*  is  in  fact  the  adjoint  of  3T  under  the  inner  product  (A 2a).  Furthermore,  the 
infinitesimal  generator  si*  of  3T*  is  the  adjoint  of  si  given  in  (2.5)  (Pazy  1983). 

The  adjoint  linearized  Navier-Stokes  equations  and  their  corresponding  evolution 
operator  form  the  basis  of  the  adjoint  input-output  system  dual  to  (2.12).  This  can 
be  obtained  in  three  steps:  (i)  derive  the  adjoint  input  and  output  operators  38*  and 
33* ;  (ii)  use  38*,  <&*  and  2T*  to  derive  the  adjoint  controllability  and  observability 
operator  and  J5?’;  and  (iii)  identify  the  adjoint  state  space  with  the  system  which 
is  associated  with  A£*c  and  jSf*. 
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The  adjoint  of  the  input  and  output  operators  38  and  ^  associated  with  the  inner 
products  (A  2 b)  and  (A  2c)  are 

r  =  (^,^)  =  (h(x;xz),  h(x;xv))  (2.25) 

and 

(®\p\  l h(x;xw)Tp\ 

@*p  =  0  =  /  0  dxdy,  (2.26) 

\®2P)  Q  \h(x;xu)TpJ 

respectively. 

The  adjoint  controllability  and  observability  operators  :  X  — *  UJ((— oo,  0])  and 
:  Y([0,  oo))  — *  X  associated  with  the  inner  products  (A  2b)  and  (A  2c)  (derived 
in  §A.2  of  Appendix  A)  are  given  by 

£?*c(-t)p(x,  0)  =  &J~(-t)p(x,  0),  (2.27a) 

/•GO 

ST0t(t)=  /  dr,  (2.21b) 

Jo 

where  t  e  X  and  p  e  X.  The  first  mapping,  J5?*,  is  from  the  adjoint  state  at  time 
t  =  0  onto  a  signal  in  HJ  at  time  —  t.  The  mapping  is  from  an  output  signal  in  X 
in  t  e  [0,  oo)  to  a  state  in  X  at  t  =  0.  In  analogy  to  the  case  of  the  forward  problem 
defined  by  (2.12),  it  can  be  seen  that  these  two  mappings  are  the  observability  and 
controllability  operators  of  the  following  state-space  system : 

—p  =  s#*p+(%*t,  (2.28a) 

e  =  STp.  (2.286) 

This  system  has  two  inputs  contained  in  the  vector  f  =  (z*,i/*)  with  t  e  X  and 
three  outputs  contained  in  the  vector  e  =  ( w* ,  u* .  g* )  e  HJ.  Comparing  the  above 
adjoint  equations  with  (2.12)  we  observe  that  the  outputs  and  inputs  have  exchanged 
places.  In  the  dual  system  (2.28),  the  adjoint  flow  field  is  forced  by  the  outputs; 
the  adjoint  problem  is  then  used  to  identify  flow  fields  yielding  the  largest  output 
response  (Dullerud  &  Paganini  1999). 

3.  Input-output  analysis 

In  this  section,  the  main  input-output  characteristics  of  our  problem  are  analysed  in 
order  to  identify  the  modes  to  retain  in  a  low-order  model.  We  introduce  the  concepts 
of  Gramians  and  balancing,  using  the  operators  defined  in  the  previous  section. 
For  a  more  detailed  presentation  of  systems  theory  please  refer  to  Kailath  (1980) 
and  Curtain  &  Zwart  (1995).  The  analysis  amounts  to  computing  the  eigenmodes  of 
three  operators:  JS?CJS?*,  0  and  jZ'cj5?*j5?*j5?£).  The  three  sets  of  modes  correspond 
to  the  flow  structures  that  are  the  most  easily  influenced  by  the  input  (controllable 
modes),  the  states  that  produce  the  largest  output  energy  (observable  modes)  and  the 
most  relevant  states  for  the  input-output  behaviour  (balanced  modes).  For  the  sake 
of  clarity,  we  will  show  numerical  results  obtained  using  only  the  first  input  tMi  and 
the  first  output  ^i,  i.e.  single-input  and  single-output  system  (SISO).  We  will  return 
to  the  multi-input  multi-output  (MIMO)  state-space  system  with  input  vector  f  and 
output  vector  y  in  §  4.  The  three  sets  of  eigenmodes  mentioned  above  can  be  computed 
numerically  for  systems  with  many  degrees  of  freedom  by  using  the  following  two 
approximations: 
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(i)  ‘The  timestepper’ :  As  mentioned  above,  solutions  of  Navier-Stokes  system  (2.12) 
in  input-output  form  are  obtained  numerically,  using  a  forward  timestepper,  which 
approximates  the  action  of  evolution  operator  .  An  adjoint  timestepper  is  used 
for  computing  solutions  of  the  associated  adjoint  system  (2.28)  and  the  action  of 
the  adjoint  evolution  operator  ST* .  The  numerical  code  employed  is  described  in 
Appendix  C.  In  the  simulations  presented,  we  have  used  768  collocation  points 
in  the  streamwise  direction  x  and  101  points  in  the  wall  normal  direction  y,  with  a 
computational  box  of  dimensions  Lx  =  1000  and  Ly  =  30  (see  figure  1).  The  discretized 
system  has  thus  m  «  105  degrees  of  freedom. 

(ii)  ‘The  snapshot  method’ :  The  controllable  and  observable  modes  introduced  next 
are  computed  using  the  snapshot  method  introduced  by  Sirovich  (1987).  Recently 
Rowley  (2005)  extended  this  method  to  obtain  balanced  modes.  The  snapshot 
technique  is  described  in  Appendix  B.  For  the  results  presented,  the  flow  structures 
are  computed  by  collecting  1600  snapshots  of  the  forward  simulation,  using  each 
input  as  initial  condition,  and  1600  snapshots  of  the  adjoint  simulation,  using  each 
output  as  initial  condition.  The  snapshots  were  taken  with  equal  spacing  in  the  time 
interval  (0,  4000). 


3.1.  Controllable  modes 

We  begin  our  input-output  analysis  by  searching  for  flow  states  that  are  most  easily 
triggered  by  a  given  input.  This  issue  is  related  to  the  concept  of  controllability,  which, 
in  general,  quantifies  the  possibility  of  steering  the  flow  between  two  arbitrary  states. 
A  state  u  is  controllable  if  it  belongs  to  the  range  of  C£c\  that  is  u  =  £Ccf(t)  exists  for 
some  f(t)f.  A  commonly  adopted  interpretation  of  controllability  is  illustrated  by  the 
following  optimal  control  problem:  what  is  the  minimum  input  energy  ||f|lu((_oo op  m 
the  time  span  t  e  (— oo,  0]  required  to  bring  the  state  (if  possible)  from  zero  to  the 
given  initial  condition  u{ x,  0)  =  m0? 

Assuming  that  u0(x)  has  a  unit  norm  and  that  it  is  controllable,  it  can  be  shownj: 
that  the  optimal  input  is  given  by 

f  =  g'^uo,  (3.1) 

where  0>  is  the  controllability  Gramian  defined  as 

/0  /»oo 

3r(-t)<M<M‘3r*(-t)dt  =  l  3r‘(t)dt.  (3.2) 

-oo  J  0 

In  the  first  equality  (2.19)  and  (2.21a)  have  been  used.  Using  (3.1)  the  minimum  input 
energy  is  given  by 

\\fWm-oo,o])=  [  uZ^-'uodxdy.  (3.3) 

JQ 

The  controllability  Gramian  3P  provides  a  means  to  rank  different  states  according 
to  how  easily  they  can  be  influenced  by  an  input.  In  particular,  the  most  easily 
influenced,  or  the  most  controllable,  flow  structures  are  the  eigenfunctions  of  SP 


f  The  system  (2.12)  is  called  exactly  controllable  if  all  states  u  e  X  can  be  reached  for  some  input. 
This  is  rarely  the  case  for  elliptic/parabolic  PDEs,  and  a  less  restrictive  condition  is  approximate 
controllability,  where  any  state  u  e  X  can  be  approximated  arbitrary  closely  by  controllable 
elements  (Curtain  &  Zwart  1995). 

j:  For  the  finite-dimensional  case  see  Lewis  &  Syrmos  (1995)  and  Dullerud  &  Paganini  (1999). 
In  the  general  case  is  well  defined  on  any  finite-dimensional  subspace  of  X. 
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Figure  3.  Instantaneous  snapshots  of  the  streamwise  disturbance  component  at 
t  =  120,  600,  1200  and  1800  triggered  by  an  impulse  in 

associated  with  the  largest  eigenvalues  of 

m  =  m-  (3.4) 

The  superscript  c  stands  for  controllable  modes.  Note  that  SP  is  a  self-adjoint 
and  positive  semi-definite  operator  whose  eigenvalues  are  real  and  positive  and 
the  eigenfunctions  mutually  orthogonal.  If  Xct  <C  1,  the  corresponding  eigenfunction 
0?  requires  very  large  energy  to  be  excited  by  the  input,  since  (A^)-1  is  proportional 
to  ||/||xu((— oo,o])-  The  mode  is  then  referred  to  as  (nearly)  uncontrollable. 

For  linear  systems  the  controllability  Gramian  corresponds  to  the  covariance  of 
the  state  response  to  an  impulse  in  time.  Therefore,  the  controllable  modes  can  be 
regarded  as  POD  modes  (Ilak  &  Rowley  2008;  Bagheri  et  al.  2008).  Traditionally, 
the  interpretation  of  these  modes  is  that  they  represent  decorrelated  energy-ranked 
flow  states.  For  example  the  first  POD  mode  <j>\  is  the  most  energetic  structure  in  the 
flow,  containing  k\/  X/i=i  '7  x  100%  of  the  total  flow  energy.  These  modes  can  be 
conveniently  obtained  by  collecting  r  snapshots  of  the  flow  at  discrete  times  t\, ...  ,tr 
and  solving  an  r  x  r  eigenvalue  problem  (Sirovich  1987). 

The  controllable  modes  can  thus  be  computed  from  the  response  of  the  flow  to  an 
impulse,  5(0): 

u(x,  tj)  =  (3.5) 

The  impulse  response  can  be  used  to  build  the  Gramian  and  compute  the  most 
controllable  modes  as  shown  in  Appendix  B.  Figure  3  shows  the  streamwise  velocity 
component  of  the  instantaneous  velocity  field  after  an  impulse  from  at  four 
different  times.  The  generation  and  convection  of  a  wavepacket  with  a  dominant 
spatial  wavenumber  and  a  propagation  speed  of  about  0.4(/oo  can  be  observed.  The 
wavepacket  grows  in  amplitude  and  size  in  the  x  direction,  until  it  reaches  the 
beginning  of  the  fringe  region  at  x  =  800.  As  it  enters  this  region,  the  disturbance  is 
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Figure  4.  The  streamwise  component  of  the  four  most  controllable  modes  0). 
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Figure  5.  (a)  The  normalized  eigenvalues  Af  (squares)  and  ?°  (circles)  associated  with  control¬ 
lable  modes  and  observable  modes  respectively.  ( b )  The  Hankel  singular  values  er;  cor¬ 
responding  to  the  balanced  modes. 


eventually  damped  by  the  fringe  forcing,  reproducing  the  effect  of  an  outflow.  The 
input-output  system  (2.12)  is  thus  asymptotically  stable. 

The  u  component  of  the  four  most  controllable  modes  0)  with  respect  to  are 
shown  in  figure  4,  while  the  corresponding  eigenvalues  A-  are  displayed  in  figure  5(a) 
with  square  symbols.  The  first  20  controllable  modes  contain  99  %  of  the  flow 
energy,  meaning  that  a  significant  part  of  the  controllable  subspace  is  spanned  by 
20  modes.  Note  that  the  flow  structures  that  are  the  most  easily  influenced  by  the 
input  are  located  downstream  in  the  domain,  where  the  energy  of  the  response  to 
forcing  is  the  largest.  In  other  words,  low  energy  is  needed  at  location  to  force 
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large  structures  downstream  owing  to  the  amplification  provided  by  the  intrinsic 
flow  dynamics.  Moreover,  the  eigenvalues  shown  in  figure  5  come  in  pairs.  The 
corresponding  velocity  fields  (see  the  first  and  the  second  modes  in  figure  4)  have 
the  same  wavepacket  structure  90°  out  of  phase.  These  modes  represent  travelling 
structures  (see  also  Rempfer  &  Fasel  1994). 

3.2.  Observable  modes 

For  a  given  sensor  it  is  important  to  determine  whether  the  relevant  flow  instabilities 
can  be  detected,  and  if  so,  to  which  accuracy.  The  flow  fields  which  can  be  most  easily 
detected  are  called  the  most  observable  modesf.  As  in  the  case  of  the  controllability 
Gramian,  the  observability  problem  can  also  be  cast  as  an  optimization  problem.  We 
wish  to  find  the  initial  conditions  that  produce  the  largest  output  energy.  The  output 
energy  generated  by  the  initial  condition  u0,  assumed  to  be  of  unit  norm,  is  given  by 

IMIy([0,oo))  =  (4? ou0i  ^oUoWd 0,oo))  =  (Ho.  25?*25?om0)x  =  j  u o  -^Xho  dx  dy,  (3.6) 

2. 

where  J  is  called  the  observability  Gramian.  Using  (2.20)  and  (2.21b)  we  obtain  the 
following  expression  for  2. : 

/••GO 

&  =  se*0se0=  3r*(t)<$*<$3r(t)dt.  (3.7) 

Jo 

The  observability  Gramian  provides  a  means  to  rank  states  according  to  their 
contribution  to  the  output.  The  most  observable  state  0"  is  given  by  the  eigenfunction 
of  the  operator  1  corresponding  to  the  largest  eigenvalue  of 

H°  =  W-  (3-8) 

The  superscript  o  stands  for  observable  modes.  Note  that  1  is  a  self-adjoint  and 
positive  semi-definite  operator  so  that  its  eigenvalues  are  real  and  positive  and 
its  eigenfunctions  mutually  orthogonal.  The  most  observable  mode  0"  contributes 
X° / XT/Lo  "A°j  x  100%  to  the  total  sensor  energy;  the  second  most  observable  mode  01,’ 
contributes  Xy  XT/=o /0j  x  100  % ;  and  so  on.  In  particular,  if  X°  <C  1,  the  corresponding 
mode  0°  does  not  make  a  contribution  to  sensor  output  and  is  called  a  (nearly) 
unobservable  mode.  Note  that  the  observable  modes  can  be  regarded  as  POD  modes 
of  the  adjoint  system. 

From  the  definition  of  Q  in  (3.7)  it  follows  that  the  observable  modes  pertaining  to  a 
given  output  can  be  determined  from  the  impulse  response  of  one  adjoint  simulation 
(see  Appendix  B).  The  results  of  this  simulation,  ,  can  then  be  used  to  build 

the  second-order  correlation  of  the  flow  field,  2r*(t)(€’‘(€2r(t),  and  thus  the  Gramian. 
The  eigenvalue  problem  (3.8)  is  solved  by  using  the  snapshot  method  as  explained 
earlier  for  the  case  of  the  controllable  modes.  Here  we  present  results  for  the  first 
output  only.  Figure  6  shows  the  instantaneous  adjoint  field  at  four  different  times 

p(x,-tj)  =  3r\tjW„  (3.9) 

after  an  impulse  from  the  first  output,  i.e.  ^(5(0).  The  triggered  wavepacket  travels 
backward  in  time  in  the  upstream  direction  with  upstream-tilted  structures.  The 
adjoint  solution  can  be  regarded  as  the  sensitivity  of  the  output  with  respect  to 
linear  perturbations  to  the  underlying  base  flow.  In  other  words,  the  flow  structures 

f  The  system  (2.12)  is  approximately  observable  if  Sfau=  0  occurs  only  when  i<=0,  i.e.  if  the 
knowledge  of  the  output  determines  the  initial  state  uniquely  (Curtain  &  Zwart  1995). 
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Figure  6.  Instantaneous  snapshots  of  the  streamwise  disturbance  component  at 
t  =  — 120,  —600,  —1200  and  —1800  of  the  adjoint  impulse  in 
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Figure  7.  The  streamwise  velocity  component  of  the  four  most  observable  modes  0". 

excited  by  and  shown  in  figure  6  are  also  the  structures  to  which  the  sensor  is 
the  most  sensitive.  In  this  context,  the  negative  time  can  be  interpreted  as  the  delay 
between  the  time  these  structures  are  present  and  the  instant  they  can  be  measured. 

The  u  component  of  the  four  most  observable  modes  0°  with  respect  to  is  shown 
in  figure  7,  while  the  associated  eigenvalues  are  reported  as  circles  in  figure  5(a).  From 
the  latter  figure,  we  observe  that  the  leading  20  modes  are  responsible  for  nearly 
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the  entire  output  energy.  The  flow  structures  in  figure  7  are  initial  conditions  that 
contribute  with  the  most  energy  to  the  sensor  output.  These  modes  are  real-valued 
functions,  and  therefore  two  of  them  are  needed  to  describe  travelling  flow  structures, 
which  explains  the  appearance  of  pairs  of  eigenvalues  in  figure  5(a).  There  are  two 
further  noteworthy  remarks : 

(i)  The  spatial  support  of  the  observable  modes  is  far  upstream,  where  the 
sensitivity  of  the  flow  is  the  largest.  Hence,  the  most  observable  structures  are 
spatially  disconnected  from  the  most  controllable  modes.  This  spatial  separation 
is  also  observed  between  the  global  eigenmodes  of  the  linearized  Navier-Stokes 
equations  and  eigenmodes  of  the  adjoint  Navier-Stokes,  where  it  is  associated  with 
streamwise  non-normality  of  the  system  (Chomaz  2005). 

(ii)  The  most  observable  structures  are  tilted  in  the  upstream  direction,  ‘leaning’ 
against  the  shear  layer,  and  are  similar  to  the  linear  optimal  disturbances  computed 
by  Akervik  et  al.  (2008).  The  optimal  disturbance  is  the  initial  condition  maximizing 
the  perturbation  energy  over  the  entire  domain  Q  at  a  fixed  time  t  =  T.  On  the  other 
hand,  observable  modes  maximize  the  time  integral  of  the  perturbation  energy  in  the 
region  defined  by  the  output  Choosing  the  sensor  location  in  correspondence  to 
the  largest  flow  response  leads  therefore  to  the  similarity  between  linear  optimals  and 
observable  modes.  As  noticed  by  Butler  &  Farrell  (1992),  the  upstream  tilting  of  the 
optimal  initial  conditions  can  be  attributed  to  the  wall  normal  non-normality  of  the 
governing  operator;  perturbations  extract  energy  from  the  mean  shear  by  transporting 
momentum  down  the  mean  velocity  gradient  (the  so-called  Orr  mechanism). 

3.3.  Balanced  modes 

So  far  we  have  identified  modes  that  characterize  either  the  response  to  forcing  or 
the  sensitivity  of  an  output.  In  this  section  we  present  the  balanced  modes  (Moore 
1981),  which  take  into  account  both  the  response  behaviour  and  the  output  sensitivity. 
Similar  to  the  previous  section,  we  wish  to  excite  the  largest  output  energy.  However, 
rather  than  identifying  dangerous  initial  conditions,  using  the  mapping  as  in  (3.6), 
we  look  directly  for  input  signals  which  produce  the  largest  output  energy  via  the 
input-output  mapping  j£?0j£?c  given  in  (2.22). 

The  output  energy  generated  by  the  past  input  f,  assumed  to  be  of  unit  norm,  is 
given  at  time  t  by 

r° 

II K II Y, [0,00))  =  ■^’^)y([0,oo))  =  {f,  3fCf) «((— oo,0])  =  /  F 2tn&t.  (3.10) 

J  —OO 

If  we  let  the  sequence  of  input  vectors  t,  with  unit  norm  represent  the  eigenfunctions 
of  i.e. 

jTjff,-  =  aftit  (3.11) 

then  the  output  energy  will  be  given  by  the  square  of  the  so-called  Hankel  singular 
values  (HSV)  ct,-.  The  most  dangerous  input  vector  f\  with  ||fi  ||u((-oo,0])  =  1  thus  results 
in  an  output  signal  which  has  been  amplified  by  erf.  Note  that  at  ^  rr2  ^  ,  so  the 

eigenmodes  of  the  input-output  map  are  ranked  according  to  how  much  the  input 
signal  is  amplified,  as  it  is  filtered  by  the  linear  system  and  the  output. 

Using  the  controllability  operator  <£ c  we  obtain  the  flow  structure  associated  with 
the  forcing  f ) : 


(3.12) 
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Notice  that  cr~1/2  is  a  convenient  normalization  factor.  The  modes  are  denoted  by  the 
superscript  oc,  which  indicates  that  these  modes  are  both  observable  and  controllable. 
The  sequence  of  functions  <j)°c  are  called  the  balanced  modes,  and  as  we  show  next, 
they  diagonalize  the  observability  Gramian.  Computing  the  output  energy  for  f,  and 
using  (3.12),  we  obtain 


fi  XT  Midi  =  >«(/,,  2?;j2C>u<(-oo,0])  =  G<c>  ^r>x  = 


=  O', 


(3.13) 


where  the  definitions  3/P  =  P£ 0£t?c,  &P*  =  and  J  =  0  are  used.  A  diagonal 

observability  Gramian  implies  that  these  modes  can  be  regarded  as  orthogonal  if 
this  Gramian  is  used  as  inner-product  weight  matrix.  With  respect  to  inner  product 
defined  in  (A  2a),  however,  these  modes  are  not  orthogonal. 

A  sequence  of  functions  f°c ,  referred  to  as  the  adjoint  balanced  modes,  which  is 
bi-orthogonal  to  according  to 

(f?'4  =  ST  <3-14) 

is  needed  to  project  our  system  on  the  basis  given  by  the  balanced  modes.  The 
derivation  is  analogous  to  <f>°c,  but  now  we  consider  instead  the  left  eigenvectors  s, 
of  the  input-output  map  J/P'J/P,  i.e. 

Jf-TTs,  =  s-,af.  (3.15) 


The  adjoint  balanced  modes  are  then  given  by 


(3.16) 


It  is  possible  to  show  by  the  same  procedure  used  in  (3.13)  that  these  modes  diagonalize 
the  controllability  Gramian: 


{tr^w 


=  07. 


(3.17) 


Furthermore,  the  diagonal  elements  are  also  equal  to  the  Hankel  singular  values. 
The  term  balancing  now  becomes  clear;  using  <j>°f  and  the  controllability  and 
observability  Gramians  become  diagonal  and  equal  to  the  HSV.  In  other  words,  the 
observability  and  controllability  properties  are  balanced.  This  is  useful  for  performing 
model  reduction,  as  it  allows  us  to  discard  the  modes  which  are  both  difficult  to 
measure  and  difficult  to  excite  by  the  inputs. 

To  compute  these  modes,  it  is  convenient  to  show  that  are  the  eigenmodes  of 
multiplying  (3.11)  with  c  yields 

=  ar-tff.  (3.18) 


The  computation  of  the  balanced  modes  and  their  associated  adjoints  can  again  be 
accomplished  using  a  timestepper  and  the  snapshot  method  described  in  Appendix  B. 
In  this  case  one  combines  the  sequence  of  snapshots  collected  from  the  solution  of 
the  forward  problem  (2.12)  with  a  sequence  of  snapshots  collected  from  the  adjoint 
system  (2.28).  In  this  way  we  can  approximate  the  eigenvalue  problem  (3.18)  to  obtain 
the  balanced  modes  (Rowley  2005).  The  u  component  of  the  first  four  balanced  modes 
<j)°c  with  respect  to  3S\  and  'Py  are  shown  in  figure  8  and  the  corresponding  adjoint 
modes  x/r°c  in  figure  9.  The  HSV  o,  are  shown  in  figure  5(h).  As  in  the  case  of  the 
observability  and  controllability  eigenvalues  /„•  and  the  singular  values  come  in 
pairs,  indicating  that  the  leading  balanced  modes  are  travelling  structures.  The  same 
observation  was  made  by  Ilak  &  Rowley  (2008)  for  channel  flow. 
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Figure  8.  The  streamwise  velocity  component  of  the  first  four  balanced  modes  0°c. 


x 


Figure  9.  The  streamwise  velocity  component  of  the  adjoint  balanced  modes  i/r°c. 

From  figures  8  and  9  we  observe  that  the  leading  balanced  modes  appear  also  as 
wavepackets,  but  they  are  somewhat  more  spatially  extended  than  the  controllable 
POD  modes  (figure  4).  Similarly,  the  adjoint  balanced  modes  have  a  larger  spatial 
support  than  the  observable  modes  (figure  7).  As  noticed  by  Ahuja  et  al.  (2007) 
and  Ilak  &  Rowley  (2008),  we  can  account  for  both  controllability  and  observability 
through  the  non-orthogonality  of  the  balanced  modes.  In  the  two  previous  sections  we 
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observed  that  for  an  input  located  upstream  and  an  output  located  downstream, 
the  associated  controllable  and  observable  modes  are  spatially  located  in  different 
parts  of  the  domain.  The  controllable  subspace  and  the  observable  subspace  are  thus 
separated  in  the  streamwise  direction.  This  is  a  consequence  of  the  convective  nature 
of  the  instabilities  arising  in  the  Blasius  flow  in  which  disturbances  grow  in  amplitude, 
as  they  are  convected  in  the  downstream  direction.  Essentially,  this  separation  implies 
that  the  distribution  of  both  the  input  and  the  output  cannot  be  captured  by  an 
orthogonal  projection  onto  the  leading  modes  of  only  one  subspace.  Conversely,  in  a 
bi-orthogonal  projection  the  adjoint  balanced  modes  account  for  the  output  sensitivity 
and  the  direct  balanced  modes  for  the  most  energetic  structures. 


4.  Model  reduction 

Since  the  disturbances  are  represented  by  an  input,  and  the  objective  consists  of 
minimizing  an  output  signal,  capturing  the  input-output  behaviour  of  the  system  - 
described  by  the  mapping  J5?0J5?C  -  is  sufficient  for  the  design  of  optimal  and  robust 
control  schemes.  The  flow  structures  that  are  neither  controllable  nor  observable 
are  redundant  for  the  input-output  behaviour.  Moreover,  the  states  that  are  nearly 
uncontrollable  and  nearly  unobservable  can  be  discarded,  since  they  have  a  very  weak 
influence  on  the  input-output  behaviour.  A  systematic  approach  of  approximating 
the  system  given  by  (2.1)  with  a  finite-dimensional  model,  which  preserves  the  main 
input-output  behaviour,  is  called  balanced  truncation  (Moore  1981).  As  we  show 
next,  balanced  truncation  amounts  to  a  projection  of  state-space  system  (2.12)  on  the 
leading  balanced  modes. 

We  now  return  to  the  MIMO  state-space  system  with  input  vector  f  and  output 
vector  y.  The  measurement  noise  acts  on  the  output  signal  and  affects  the  perturbation 
dynamics  only  in  the  closed-loop  system  and  is  hence  not  included  in  the  analysis. 

4.1.  Galerkin  projection 

Any  flow  field  can  be  approximated  as  a  linear  combination  of  the  leading  r  balanced 
modes : 

r 

ur(x,t)  =  j2^mojc(x),  (4.i) 

j  i 

where  c/j  =  {u,  is  the  expansion  coefficient.  Inserting  the  above  expansion 

into  (2.12)  and  taking  the  inner  product  with  the  adjoint  balanced  modes  x/r°c , 


the  following  r -dimensional  state-space  form  is  obtained: 

q  =  Aq  +  B\W  +  B2u,  (4.2  a) 

v  =  Ciq+lu,  (4.2  b) 

z  =  C2q  +  ocg  .  (4.2 c) 

This  system  is  referred  to  as  the  ROM.  The  column  vector  q  contains  qj,  and  the 
entries  of  the  matrix  A,  column  vector  B\  and  row  vector  Cj  are 

A,- (4.3a) 
B1J  =  (^C,^1>X,  (4.36) 

Cij  =  Vi  <t>°f  (4.3c) 


282 


S.  Bagheri,  L.  Brandt  and  D.  S.  Henningson 


Figure  10.  (a)  The  Flankel  singular  values  (black  symbols)  are  compared  to  the  diagonal 
entries  of  the  controllability  (red)  and  observability  (blue)  Gramians  associated  with  the 
balanced  reduced-order  system.  ( b )  The  model  reduction  error.  The  upper  and  lower 
theoretical  bounds  are  depicted  with  grey  lines,  and  the  actual  model  reduction  error  is  shown 
with  black  symbols. 

for  2,  j  =  1, . . . ,  r.  The  components  of  the  vectors  C2  and  B2  are  obtained  in  the  same 
manner  as  those  of  B\  and  C\.  The  evolution  operator  associated  with  (4.2)  is 

T(t)  =  eA'  =  YJ-^--  (4.4) 

j= o  J ' 

Notice  that  the  balanced  modes  are  computed  accounting  for  all  the  inputs  (except 
the  measurement  noise)  and  outputs,  and  the  Galerkin  projection  (4.2)  is  performed 
only  once.  The  projection  of  si  on  the  balanced  modes  can  be  approximated  by  the 
finite-difference  method,  using  the  timestepper  and  (2.5).  For  the  results  presented,  8t 
was  chosen  to  be  10  4  after  a  convergence  study. 

To  validate  the  properties  of  the  snapshot-based  balanced  truncation,  we 
constructed  the  reduced  model  (4.2)  and  computed  its  controllability  and  observability 
Gramians.  The  projected  system  is  internally  balanced  only  if  its  Gramians  are 
diagonal  and  equal  to  the  HSV.  We  found  that  the  first  70  x  70  elements  of  both 
Gramians  were  diagonal.  In  figure  10  we  compare  the  leading  100  diagonal  elements 
with  the  HSV.  The  first  70  modes  are  observed  to  be  bi-orthogonal  to  each  other 
down  to  numerical  accuracy.  However,  for  higher  modes,  as  the  numerical  round-off 
errors  increase,  the  bi-orthogonality  is  gradually  lost,  and  off-diagonal  elements  are 
observed  in  both  Gramians.  By  increasing  the  numerical  resolution  and  the  number 
of  snapshots  it  is  possible  to  increase  the  number  of  balanced  modes.  However,  -  as 
noticed  by  (Moore  1981)  -  the  ratio  cti/ct,  serves  as  a  condition  number  for  fif,  and 
therefore  the  balanced  modes  corresponding  to  very  small  HSV  can  be  ill  conditioned 
independent  of  the  numerical  approximations. 

4.1.1.  Performance  of  the  ROM 

In  this  section  the  input-output  behaviour  of  the  ROM  (4.2)  is  compared  to  the  full 
Navier-Stokes  system  (2.12).  We  begin  by  comparing  the  impulse  response  from  all 
inputs  to  all  outputs.  In  figure  11  three  signals  ddi  — >  ^i,  — >  (€2  and  Sd2  — >  are 

shown  with  black  lines.  The  response  of  (&2  to  forcing  in  &2  is  zero,  since  disturbances 
travelling  upstream  are  quickly  damped.  These  impulse  responses  were  obtained  by 
using  the  timestepper  with  ~105  degrees  of  freedom.  The  impulse  responses  of  the 
ROM  (4.2)  with  r  =  50  given  by  y(t)  =  CeAl  B  are  shown  with  red  dashed  lines.  We 
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Figure  11.  The  impulse  response  from  38\  — >  2  (a),  38 2  — >  ^1  (b)  and  3fl\  — *  'g’l  (c).  The 

black  solid  line  represents  direct  numerical  simulations  with  105  degrees  of  freedom  and  the 
red  dashed  line  the  balanced  reduced-model  with  50  degrees  of  freedom. 

observe  that  reduced  model  registers  the  same  signal  as  the  full  model  from  all 
inputs  to  all  outputs.  The  wavepacket  triggered  by  the  impulse  of  38\  reaches  the 
first  sensor  after  600  time  units  and  the  second  sensor  (€\  after  1500  time  units. 
The  wavepacket  triggered  from  the  actuator  reaches  the  second  sensor  after  600  time 
units. 

The  frequency  response  of  the  full  system  and  the  ROM  are  compared  next. 
The  frequency  response  is  related  to  the  Laplace  transform  of  the  impulse  response 
(08  — )  (see,  e.g.  Skogestad  &  Postlethwaite  2005),  which  in  our  case  results  in  the 
2x3  transfer  function  matrix  (TFM): 

G(s)  =  <€(s  -  (4.5) 

with  ss€.  The  element  G,j  contains  the  response  from  38 j  — >  (€i.  The  TFM  of  size 
2  x  3  of  the  reduced  model  is  similarly  defined  as 

Gr(s)  =  C(sl  —  A)~lB  (4.6) 

with  1  as  identity  matrix  of  size  r. 

Due  to  the  linear  nature  of  the  equations,  a  sinusoidal  input  signal  e10Jf  with 
constant  frequency  co  will  generate  an  output  with  the  same  frequency  but  with  a 
phase  shift  Arg[G(ia>)}  and  a  different  amplitude  |G(iw)|.  The  frequency  response 
gain  is  usually  measured  by  the  largest  singular  value  of  the  TFM  (Skogestad  & 
Postlethwaite  2005).  For  the  full  model  we  do  not  have  an  explicit  expression  of  the 
TFM.  Therefore,  we  make  use  of  our  timestepper  and  apply  a  sinusoidal  signal  with 
a  constant  frequency  co  in  each  input  and  extract  the  periodic  signal  from  the  outputs 
once  the  initial  transients  have  died  out.  Note  that  computing  the  frequency  response 
with  the  timestepper  in  this  way  does  not  take  into  account  the  interaction  of  input 
signals,  since  only  one  input  is  active  at  a  time. 


284 


S.  Bagheri,  L.  Brandt  and  D.  S.  Henningson 


Figure  12.  The  envelope  of  the  MIMO  transfer  function  matrix  G(io>)  from  all  inputs  to  all 
outputs  computed  using  the  timestepper  is  shown  with  red  symbols.  The  largest  response  is  for 
a>  =  0.06  with  a  peak  value  of  144.6.  For  w  e  [0,  0.03]  the  frequency  response  from  the  actuator 
to  objective  function  — ►  @i)  dominates.  The  frequencies  with  the  largest  gain  are  obtained 

from  disturbances  to  objective  function  (dSi  —>  %{)  in  the  range  co  e  [0.03,0.07],  and  finally 
for  higher  frequencies  the  response  from  disturbances  to  measurement  sensor  (^i  — >  ‘gh)  are 
amplified  the  most.  The  frequencies  in  the  grey  domain  are  amplified.  Also  shown  are  the 
frequency  response  (the  envelope)  of  the  reduced  model  TFM  Gr(ia>)  with  ranks  2  (green),  50 
(blue)  and  100  (black). 


In  figure  12  the  envelope  of  the  TFM  amplitudes  -  the  largest  amplification  of  all 
the  frequency  responses  from  tffii  — >  #1,  — >  (&1  and  — ►  ^1  at  each  co  -  for  the 

full  model  of  order  105  is  shown  with  red  filled  circles.  In  the  same  figure  the  TFM 
amplitudes  of  ROMs  of  order  r  =  2,  50  and  100  are  shown.  We  observe  that  the  ROM 
of  order  2  captures  the  most  important  aspect  of  the  input-output  behaviour,  which  is 
the  response  of  the  most  dangerous  frequency,  i.e.  the  peak  response  of  the  full  model. 
The  model  with  50  modes  is  able  to  estimate  the  gains  of  all  the  amplified  frequencies 
but  fails  to  capture  the  damped  low  and  high  frequencies.  Adding  50  additional 
modes  results  in  a  model  that  preserves  the  input-output  behaviour  correctly  for 
all  frequencies.  Note  that  there  are  no  isolated  eigenvalues  in  the  spectrum  of  the 
spatially  developing  Blasius  flow  (Ehrenstein  &  Gallaire  2005;  Akervik  et  al  2008), 
and  therefore  the  frequency  response  is  rather  smooth  with  no  peaks.  Low-pass  filters 
of  this  form  cannot  be  represented  with  only  a  few  degrees  of  freedom. 

Finally,  the  model  reduction  error  is  computed  and  compared  to  the  theoretical 
bounds  given  by  the  Hankel  singular  values.  An  attractive  feature  of  balanced 
truncation  is  the  existence  of  error  bounds  (which  are  obtained  a  priori  to  Galerkin 
projection): 

n 

CT,.+1  <  II G  —  Gr||oo  <  2  ^2  17 J-  (4.7) 

j=r+ 1 

The  infinity  norm  of  the  transfer  function  equals  the  peak  value  of  the  frequency 
response  gain  based  on  the  largest  singular  value,  i.e. 

||G(j)||oo  =  maxffi(G(ifl))).  (4.8) 

CO 

Estimating  the  model  reduction  error  (4.7)  amounts  to  the  calculation  of  the  difference 
of  the  peak  values  of  frequency  response  of  the  reduced-order  and  the  Navier-Stokes 
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Figure  13.  The  closed-loop  system.  The  plant  represents  the  input-output  system  given 
by  (2.12)  subject  to  external  disturbances  w.  The  controller  of  the  order  50  forces  the 
Navier-Stokes  equations  with  the  input  signal  u  based  on  the  noisy  measurements  v,  so 
that  the  effect  of  w  on  the  output  signal  z  is  minimized. 

systems.  For  the  latter  system,  we  use  the  peak  value  of  the  amplitude  envelope  as 
shown  in  figure  12  instead  of  the  largest  singular  value.  The  error  norm  for  the 
balanced  truncation  model  is  shown  in  figure  10(h)  with  black  symbols.  The  error 
norm  is  close  to  the  lower  bound  given  by  the  HSV  for  the  first  50  modes.  The  peak 
value  for  the  Navier-Stokes  system  is  144.6  which  is  gradually  approached  by  the 
ROM,  until  it  saturates  at  a  peak  value  of  144.5  due  to  numerical  round-off  errors. 
Note  that  the  error  is  somewhat  lower  than  the  theoretical  bounds  for  the  reduced 
systems  of  orders  2  and  4.  This  is  because  the  frequency  response  of  the  full  system 
is  obtained  numerically  using  our  timestepper  and  because  HGHo,  is  based  on  the 
maximum  of  the  envelope  of  the  TFM  instead  of  its  largest  singular  value. 

A  thorough  comparison  between  the  ROMs  obtained  with  POD  modes  and 
balanced  modes  can  be  found  in  Ilak  &  Rowley  (2008)  for  the  case  of  channel 
flow  and  Bagheri  et  al.  (2008)  for  the  linear  Ginzburg-Landau  equation.  The  latter 
work  also  includes  global  eigenmodes  of  the  linearized  operator  for  comparison. 

5.  Feedback  control 

We  will  now  develop  a  reduced-order  feedback  controller,  which  will  have  the  same 
dimension  as  the  ROM  (e.g.  r  =  50).  The  closed-loop  behaviour  and  the  objective  func¬ 
tion  z  will  be  investigated  and  compared  to  the  uncontrolled  flat-plate  boundary  layer. 

5.1.  3^2 -framework 

The  main  idea  of  linear  feedback  control  is  shown  in  figure  13.  As  stated  in  the 
introduction  our  objective  is  to  find  a  control  signal  u(t),  such  that  in  the  presence 
of  disturbances  w(t),g(t)  the  perturbation  energy  represented  by  the  state  variable 
u( x,  t)  is  minimized  downstream  at  the  location  defined  by  the  sensor  This  is  the 
control  problem. 

In  the  previous  section  we  showed  that  our  reduced  model  (4.2)  is  able  to  capture 
the  input-output  behaviour  of  the  Navier-Stokes  system  (2.12).  During  the  control 
design  process  we  can  assume  that  the  reduced  model  is  the  plant  that  we  wish  to 
control.  Once  we  have  determined  the  control  law  for  this  approximating  model,  we 
will  apply  it  to  the  full  Navier-Stokes  system.  Please  refer  to  Anderson  &  Moore 
(1990),  Zhou  et  al.  (2002)  and  Bagheri  et  al.  (2008)  for  further  details  of  the  Jf2 
control  algorithm,  as  we  will  only  outline  the  main  steps. 
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Following  the  notation  introduced  for  the  reduced  model  (4.2),  the  objective 
function  (2.11)  becomes 

/»00 

IIz!Iy([0,oo))  =  /  qTC{Ciq +12UTU  dr.  (5.1) 

Jo 


The  determination  of  the  control  signal  is  based  only  on  the  measurements  from 
the  sensor  f£2.  However,  for  linear  systems  -  due  to  the  separation  principle  (Zhou 
et  al.  2002)  -  the  feedback  control  law  can  be  determined  assuming  that  the  complete 
velocity  field  is  known.  The  forcing  needed  to  reproduce  the  flow  only  from  wall 
measurements  can  be  computed  independently.  Hence,  the  control  design  of  the  Jf2 
control  is  performed  in  the  following  three  steps: 

(i)  Compute  the  control  feedback  gain  K  by  solving  a  Riccati  equation  (see 
Appendix  D),  so  that  the  control  signal  is  of  feedback  type,  i.e. 


u(t)  =  Kq(t). 


(5.2) 


This  leads  to  a  new  system  (compared  to  (4.2)) : 


q  —  (A  T  B2K)q  -f*  B\W,  (5.3u) 

z  =  Ciq.  (53b) 


It  is  expected  that  the  above  perturbated  operator  A  +  B2K  has  dynamics  that  result 
in  a  smaller  amplitude  of  the  output  signal  z  than  for  the  unperturbated  operator  A 
in  (4.2). 

(ii)  Compute  the  estimation  feedback  gain  L  also  by  solving  a  Riccati  equation  (see 
Appendix  D),  so  that  the  observer 

q  —  (A  T  LC2)q  T  Lv  (5.4u) 

is  asymptotically  stable,  i.e.  \\q  —  q\\  — >■  0  as  t  — >  oo.  This  implies  that  the  estimated 
state  q  based  on  the  measurements  v  approaches  the  true  state  q  exponentially  fast. 

(iii)  The  compensator  (controller  in  figure  13)  is  finally  obtained  as 

q  =  (A  T  B2K  T  LC2)q  —  Tv,  (5.5 a) 

u  =  Kq.  (5.5  b) 


Given  the  measurements  signal  v  from  the  physical  flow,  the  reduced-order  controller 
provides  an  optimal  control  signal  u  proportional  to  the  estimated  flow  q. 

To  apply  feedback  control  in  the  numerical  simulations,  an  augmented  state-space 
system  with  state  (u,q)T  is  considered:  its  dynamics  are  given  by  (2.12)  and  (5.5), 
inputs  ( w ,  g )  and  with  the  single  output  z: 


(u\  _  (  si  @2K  \  (u\  ( °\(w\ 

\ A  +  B2K  +  LC2)  +  ^  0  -L)\g)' 

z  =  ^xU. 


(5.6  a) 
(5.6  b) 


This  system  is  referred  to  as  the  closed-loop  system.  Note  that  the  feedback  gain  K 
and  estimation  gain  L  have  the  dimension  of  the  reduced  model,  resulting  in  a  fast 
online  controller. 

The  spatio-temporal  evolution  of  the  perturbations  governed  by  the  closed-loop 
system  is  obtained  by  solving  (5.6)  numerically  using  the  timestepper  described  in 
Appendix  C  and  the  small  reduced  system  in  (5.5)  simultaneously.  The  latter  system 
is  solved  using  a  standard  Crank-Nicholson  scheme. 


Input-output  analysis,  model  reduction  and  control  of  flat-plate  boundary  layer  287 


(c) 


100 

o 

-100 


(d) 


0  1000  2000  3000  4000  5000  6000  7000  8000  9000  10  000 


t 

Figure  14.  Input  and  output  signals  of  the  closed-loop  system.  The  (a)  random  forcing  w, 
( b )  measurements  signal  v ,  (c)  control  signal  u  and  ( d )  objective  function  z  are  shown.  The 
cheap  controller  is  active  between  t  e  [2500,  7500],  represented  by  the  grey  area. 


5.2.  Performance  of  closed-loop  system 

We  will  now  investigate  the  performance  of  the  closed-loop  system  (5.6).  In  particular, 
the  output  z  of  the  closed-loop  -  with  optimal  control  signal  u  -  and  of  the  linearized 
Navier-Stokes  equations  without  control  are  considered  in  the  case  of  stochastic  and 
harmonic  forcings  in  w. 

Three  controllers  are  investigated:  (i)  cheap  control/low  noise  contamination  with 
/  =  0.1  and  a  =  0.1;  (ii)  expensive  control/high  noise  contamination  with  /  =  10  and 
a  =  10;  and  (iii)  an  intermediate  case  with  /  =  2  and  a  =  0.1. 

Note  that  the  purpose  of  the  measurement  noise  g  is  to  account  for  uncertainties  in 
the  sensor  measurements  during  the  control  design.  When  evaluating  the  closed-loop 
performance  -  solving  the  controlled  Navier-Stokes  equations  -  the  system  is  only 
forced  with  w  and  not  with  g. 

The  performance  of  the  control  in  case  (i)  is  examined  first.  In  figure  14  the 
input  and  output  signals  are  shown.  The  grey  region  indicates  the  time  at  which  the 
control  is  active.  As  disturbance  signal  w{t)  we  choose  white  noise;  the  corresponding 
response  of  the  sensor  v(t)  in  figure  14(h)  confirms  the  amplification  and  filtering 
of  the  signal  as  it  traverses  the  unstable  domain.  The  disturbances  reach  the  second 
sensor  (figure  14 d)  after  about  1500  time  units,  where  they  are  amplified  by  one  order 
of  magnitude.  The  control  is  activated  at  time  t  =  2500;  the  actuator  immediately 
begins  to  force  the  system  with  a  control  signal  (figure  14c)  based  on  the  output  v, 
and  after  a  delay  of  another  1500  time  units,  the  stabilizing  effect  of  the  control  signal 
on  the  output  z  is  clear.  When  the  control  is  deactivated  (at  t  =  7500)  the  disturbances 
start  to  grow  again. 

The  wall  normal  maximum  of  the  root  mean  square  (r.m.s.)  values  of  the  streamwise 
velocity  component  in  cases  with  and  without  control  are  shown  in  figure  15.  The 
r.m.s.  value  grows  exponentially  downstream  in  the  uncontrolled  case  until  the  fringe 
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Figure  15.  The  r.m.s.  values  of  the  uncontrolled  system  (red  line),  cheap  controller  (solid 
black),  intermediate  controller  (dashed-dotted  line)  and  expensive  controller  (dashed  line). 
The  grey  bar  represent  the  size  (defined  as  99  %  of  the  spatial  support)  and  the  location  of 
the  two  inputs,  whereas  the  red  bars  correspond  to  the  two  outputs. 


Figure  16.  Comparison  of  the  frequency  response  from  disturbances  to  objective  function 
— >  'g’l)  of  open-loop  (red)  and  three  closed-loop  systems.  The  cheap,  intermediate  and 
expensive  controllers  are  represented  by  the  solid  black,  dashed-dotted  and  dashed  lines 
respectively.  The  infinity  norm  of  the  open  loop  is  ||G|oo  is  140.7,  whereas  for  the  closed-loop 
systems  ||Gc||oo  it  is  6.4  (cheap),  9.4  (intermediate)  and  101.9  (expensive). 


region  at  x  =  800.  The  r.m.s.  of  the  controlled  perturbation  grows  only  until  it  reaches 
the  actuator  position  at  which  it  immediately  begins  to  decay.  At  the  location  of  the 
objective  function  ^  (x  =  750),  the  amplitude  of  the  perturbations  is  one  order  of 
magnitude  smaller  than  in  the  uncontrolled  case  for  the  cheapest  controller. 

The  r.m.s.  values  in  the  case  of  the  expensive  (case  ii)  and  intermediate  (case  iii) 
controls  are  shown  with  dashed  and  dashed-dotted  lines  respectively.  The  expen¬ 
sive  control  is  very  conservative,  as  the  measurement  signals  are  highly  corrupted 
and  the  control  effort  limited;  it  results  only  in  a  small  damping  of  the  disturbances. 
The  intermediate  controller  (case  iii)  is  more  cautious  in  reducing  the  perturbation 
energy  just  downstream  of  the  actuator  when  compared  to  the  cheap  controller. 
However,  it  is  interesting  to  note  that  at  the  location  at  which  the  objective  function 
is  measured,  the  disturbance  amplitude  has  decreased  nearly  as  much  as  with  the 
cheap  controller,  although  the  total  perturbation  energy  is  much  larger  over  the  entire 
domain. 

In  figure  16  the  frequency  response  from  w->z  of  the  uncontrolled  Navier-Stokes 
equations  (2.12)  (shown  in  red)  is  compared  to  that  pertaining  to  the  three  controllers 
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under  consideration.  The  solid  black  line  corresponds  to  cheap  control,  dashed- 
dotted  line  to  intermediate  control  and  dashed  line  to  expensive  control.  The  first  two 
controllers  suppress  the  most  dangerous  frequencies  close  to  co  =  0.6  significantly.  Note 
that  compared  to  the  uncontrolled  model,  the  highly  damped  frequencies  &>>0.11 
have  larger  gain  in  amplitudes.  This  behaviour  is  often  observed  in  closed-loop 
physical  systems  and  is  related  to  the  ‘waterbed’  effect;  i.e.  when  certain  frequencies 
are  suppressed,  the  response  at  other  frequencies  is  amplified. 


6.  Conclusions 

Model-based  feedback  control  of  the  instabilities  arising  in  a  spatially  inhomo¬ 
geneous  boundary  layer  flow  is  studied.  To  build  a  reduced-order  model  of  the 
problem,  where  the  application  of  standard  tools  from  control  theory  become 
computationally  feasible  also  for  fluid  flow  systems,  the  main  features  of  the  flow 
behaviour  are  investigated  in  an  input-output  framework.  The  observable,  controllable 
and  balanced  modes  of  the  system  have  been  identified.  The  location  and  structure  of 
these  modes  reflect  the  location  of  sensors/actuators  and  the  perturbation  dynamics; 
i.e  the  observable  modes  are  located  upstream,  where  the  sensitivity  to  initial 
conditions  is  the  largest.  The  controllable  modes,  conversely,  are  located  downstream, 
where  the  response  to  the  forcing  is  the  largest.  The  analysis  presented  here  can  be 
closely  related  to  stability  analysis,  using  global  modes  and  optimal  disturbances, 
except  that  inputs  and  outputs  are  taken  into  account.  The  quantity  one  wishes  to 
optimize  for  is  now  defined  by  a  sensor  output,  while  perturbations  are  introduced  by 
the  inputs  considered  in  the  model.  Furthermore,  in  view  of  the  control  application, 
the  formulation  of  the  control  objective  function  as  an  output  is  particularly 
attractive  in  this  input-output  setting,  since  this  behaviour  is  well  captured  by  the 
ROM. 

Model  reduction  is  achieved  by  projecting  the  governing  equations  on  the  leading 
balanced  modes  of  the  system.  We  show  that  the  input-output  behaviour  of  the 
flat-plate  boundary  layer  can  be  captured  accurately  with  an  ROM  based  on  these 
modes.  Finally,  the  model  is  used  to  apply  feedback  control  based  on  measurements 
from  one  upstream  sensor  and  an  actuator  further  downstream.  The  perturbations 
growth  could  be  reduced  efficiently  using  the  optimal  feedback  controller. 

It  is  also  important  to  note  that  the  approach  followed  here  requires  only 
the  use  of  a  timestepper,  a  numerical  code  solving  the  Navier-Stokes  equations, 
and  avoids  the  use  of  the  large  matrices  defining  the  operators  governing  the 
input-output  behaviour.  In  addition,  the  present  formulation  accounts  naturally  for 
localized  sensors  and  actuators,  and  therefore,  it  can  be  directly  applied  to  different 
flow  configurations.  We  are  currently  extending  the  analysis  to  three-dimensional 
disturbances  and  also  incorporating  more  realistic  actuators  (blowing/suction)  and 
sensors  (wall  measurements).  These  computations  are  now  feasible  and  will  take  us  a 
step  closer  to  using  the  controller  in  actual  experiments. 
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Appendix  A.  Derivation  of  adjoint  operators 

A.l.  The  adjoint  operators  s/*,  tffl*  and  c€* 

For  a  bounded  linear  operator  between  two  Hilbert  spaces,  J5?  :  Xi  — >  X2,  the  adjoint 
operator  satisfies 


(g’q,  p)Xl  =  (q,  p)^  for  all  (Al) 

The  derivations  make  use  of  the  following  definitions  of  inner  products : 


(■ u,p)x=  /  u{x)T  p(x)dx  dy, 
(f,  g)v  =fTg , 

(z,y)Y  =zTy , 


Vh,  p  e  X, 
W,  g  e  UJ, 
Vz,  /e¥, 


>  (A  2) 


it  g)u((-oo,o]) 
(z-  y)Y([o,oo)j 


r  t 

=  /  zT  y  dt, 


Vf  g  e  TU(( — oo,  0]), 
Vz,  ye  Y([0,oo)). 


jo  t 

Note  that  the  kinetic  energy  of  a  perturbation  u  at  time  t  is  measured  by  ||h||x  = 
( u ,  m)x-  We  begin  by  deriving  the  adjoint  operator  of  :  UJ  — >  X,  using  the  identity 

m)x  =  {f,  (A  3) 

The  left-hand  side  is  equivalent  to 


{(Mf)ru  dx  dy  =  f 


u  dx  dy  =  (f, 


u  dx  dy)iu; 


(A4) 


using  (A3)  we  identify  :  X  — >  UJ 


u 


r u  dx  dy. 


(A  5) 


The  adjoint  of  the  output  operator  ^  :  X  — >  X  can  be  derived  in  an  analogous 
manner  by  using  the  identity 


(^k,  y)x  =  {u,%*y}x-  (A 6) 

The  left-hand  side  can  be  written  as 

(j$u)T y  =  f  uT%Ty  =  (u,$>Ty)x,  (A 7) 

J  a 

where  ^  is  the  integrand  in  (2.17).  We  can  now  identify  the  adjoint  output  operator 

r:Y^Xas 

<ty  =  %Ty.  (A  8) 

The  evolution  operator  2T  :  X  — *  X  has  been  defined  in  (2.4).  The  adjoint  of  3~ 
satisfies 

(Fu,  p)x  =  {u,dT"p)x.  (A9) 

We  begin  with  taking  the  inner  product  of  p  and  a  with  the  Navier-Stokes 
equations  (2.1a)  and  (2.1/?),  respectively.  By  integrating  over  the  time  domain  and 
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applying  integration  by  parts  we  obtain 


°  =  /  /  1  P 

/o  j Q 


rp  I  dll  T  T  1  9 

T  1  —  -  (-( V  •  V)  -  (VUT)T+Re-lS/2  +  a(x))u  -  Ip 
at 


0  J  Q 


—u 


+<t(V  *u) Jdx  d y  d t 

+  ((U  ■  V)  -  (VU)T  +Re~LV2  +  a(x))p  +  Vct  )  -  p(V  ■  p)  j  dx  dy  dr 


T  I  dfl  ,  (<TJ  IX7TT\T  i  D^,-lY72 

dt 


+  /  [B.C.]^  df  +  /  u(t)T  p(t)dx  dy  —  /  M(0)rp(0)  d.r  dy 


(A  10) 


By  requiring  the  first  two  terms  to  be  zero  we  obtain  the  adjoint  Navier-Stokes 
equations  with  boundary  conditions.  They  will  be  considered  after  the  boundary  terms 
in  time  denoted  by  3  in  (A  10).  We  thus  require  that 

f  u(t)T  p(t)dx  dy  =  f  u(0)T  p(0)  dx  dy.  (All) 

The  left-hand  side  can  be  rewritten  as 

[  {^~{t)u(0))T p(t)dx  dy  =  (&~(t)u(0),  p(t))jc 
J  £2 

=  (H(0),^"*(r)p(r))x  =  [  u(d)r 2r\t)p(t) , 

J  Q 

where  we  can  identify  the  action  of  the  adjoint  evolution  operator  2T*  :  X  —>  X  as 

*~(t)p(t)  =  p(  0).  (A  12) 

Now  we  proceed  with  deriving  the  adjoint  equations  associated  with  ST* .  The  spatial 
boundary  terms  given  by  the  second  term  in  (A  10)  are 

[  [B.C.]j2  dr  =  [  cru  +  u p  +  UuT p  —  Re~l pT—  +  Re ^u7  — 

Jo  Jo  [  dx  3vJ0 

.  1  Ly 


+ 


crv  +  v* p  +  VuTp  —  Re  lpT—+Re  luT  — 

dy  dy 


dt  =0. 


J  o 


If  boundary  conditions  (2.2)  on  u  are  used  and  if  we  demand  that  p  =  (u* ,  v*),  cr* 
and  that  p  satisfy 


(cr,  p)( 0,  y)  =  (ct,  p)(Lx,  y),  (A  13a) 

p(0,y)  =  p(Lx,y),  (A  13  b) 

p(x,  0)  =  p(x,  Ly)  =  0,  (A  13c) 

we  observe  that  the  boundary  terms  vanish. 

Finally,  the  first  term  in  (A  10)  defines  the  adjoint  Navier-Stokes  equations  if  we 
demand  that  p  satisfy  (2.23).  Together  with  boundary  conditions  (A  13b)  and  (A  13c), 
(2.23)  determine  the  behaviour  of  adjoint  flow  field  p. 

A.2.  The  adjoint  operators  T£*c  and 

The  adjoint  of  the  controllability  operator  c  :  HJ(( — oo,  0])  — »  X  is  derived  using  the 
identity 


(=S? cf,  «o)x  =  (f,  =2'*«o)lJ((-oo,0])- 


(A  14) 
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We  expand  the  left-hand  side: 


(=S? c/,  «o)x  —  [  [ 
J  Q  J  —  o 

-L 


t)ddf(t))Tu0  dr  dx  dy 


=  /  dr 

J  —00 

=  (/(t),^*^*(-t)Mo)lU,(-oo,0])- 

In  the  first  equality  the  definitions  of  £%*  and  2T*  from  (A  5)  and  (A 9)  were  used.  We 
can  now  identify  :  X  — >  U((— oo,  0]) 

j5f*(-r)Ho  =  @*3r*(-t)uo.  (A  15) 

In  a  similar  fashion  the  adjoint  of  the  observability  operator  :  X  — >  Y([0,  oo)) 
is  defined  by 

y)v([o,oo))  =  (h,  (A  16) 

Expanding  the  left-hand  side  results  in 

/•GO 

{seQu,  y)v([0,OT))  =  /  (^(t)M(f))rydr 

JO 


fT  (t)u(t))T  y  dx  dy  dr 


/0 


u  (3T* (t)^* y(t))  dx  dy  dr 


/0 


=  <  U , 


rrydr 


where  is  the  integrand  in  (2.17).  We  can  identify  the  adjoint  observability  operator 
ST0  :X([0,oo))^Xas 


/•OO 

Kvif)=  /  2T\t)ry{t)dt. 
Jo 


(A  17) 


Appendix  B.  The  snapshot  method 

We  will  show  how  to  approximate  the  eigenvalue  problems  (3.8),  (3.4)  and  (3.18)  in 
order  to  compute  the  observable,  controllable  and  balanced  modes. 

B.l.  Approximate  Gramians 
We  begin  with  considering  the  eigenvalue  problem 

P+C,=W'  (Bl) 

where 

/•GO 

&#=  /  (B  2) 

Jo 

The  first  step  is  to  rewrite  the  action  of  the  controllability  Gramian  SP  on  in  terms 
of  impulse  responses  of  the  state.  Recall  that  the  flow  field  triggered  by  an  impulse 
5(0)  applied  to  the  input  is  given  by  ST Let  us  define  the  vector  u  containing 
the  impulse  responses  of  all  inputs  (except  the  measurement  noise)  as  columns,  i.e. 

ft  =  2T 3ft  =  (^dx),  .TdSjix))  =  («i(jc,  r),  u2{x ,  r)).  (B  3) 
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Moreover,  from  the  expression  of  and  2T*  given  by  (A  5)  and  (A  9)  respectively, 
we  can  rewrite  the  action  of  tffl*  2T*  on  as 


—  /  t art&>\T a.c 


(ray  </>■  dx  dy 


«r0-  dx  dy. 


(B4) 


The  controllability  Gramian  becomes 


r 


uT dx  dy  I  df. 


(B  5) 


This  expression  is  approximated  by  discretization  in  space  and  time.  Suppose  that 
x  =  (jti, . . . ,  xn/2)  is  a  grid  in  Q  with  n/2  points.  We  construct  an  n  x  2  matrix  by 
evaluating  u  at  the  grid  points,  i.e. 


/  Ui(X\,t)  U2(x  1 ,  t ) 

•  5 

\Ui(xn/2,t)  U2(xn/2,t) 


(B  6) 


The  size  of  this  matrix  is  n  x  2  because  it  has  two  velocity  components,  e.g.  uflxj,  t)  = 
(«!_!,  (xj,  t),  uy2(x  j,  t))T .  Similarly  we  define  0,  as  the  following  nxl  column  vector: 


0,  =  (^(xi),...,^(xB/2))r.  (B7) 

The  integral  in  £2  in  (B  5)  can  now  be  approximated  by 


/  u7^  dx  dy  «  uTW(j)i,  (B8) 

Jn 

where  the  n  x  n  positive-definite  matrix  W  contains  the  spatial  integration  weights 
8X  .  The  quadrature  weights  8Xj  depend  on  the  chosen  quadrature  rule.  For  instance 
in  our  case,  8X.  consist  of  the  Chebyshev  integral  weight  functions  (Hanifi,  Schmid 
&  Henningson  1996)  in  the  wall  normal  direction  and  a  trapezoidal  rule  in  the 
streamwise  direction. 

The  expression  given  by  (B  5)  becomes 


W4>i, 


(B  9) 


where  we  recognize  the  term  in  the  parenthesis  as  the  state-covariance  matrix.  If 
the  flow  fields  are  given  as  snapshots  at  discrete  times  t\, . . .  ,tm,  we  can  further 
approximate  (B9)  with 


xxTwd>i. 


(BIO) 


The  n  x  2m  matrix  X  contains  u(tj)  in  column  /,  multiplied  with  the  square  root  of 
the  quadrature  coefficients  8t  in  time,  i.e. 


A  =  (uff^/S^, . . . ,  u(tm)s/^J,  (B  11) 


where  each  column  of  X  is  referred  to  as  a  snapshot. 

The  eigenvalue  problem  given  by  (B  1)  can  now  be  approximated  with  the  following 
n  x  n  eigenvalue  problem: 

XXTWfi  =A^i, 


i  =  1, ...  ,n. 


(B  12) 
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It  is  prohibitively  expensive  to  diagonalize  the  matrix  XXTW  when  n  ^  105.  In 
the  method  of  snapshots  (Sirovich  1987),  the  modes  0,-  can  be  approximated  by 
diagonalizing  the  2m  x  2m  matrix  XTWX  instead.  This  is  efficient  when  the  product 
of  the  number  of  snapshots  and  the  number  of  inputs  is  much  smaller  than  the 
number  of  spatial  grid  points. 

In  the  method  of  snapshots  the  modes  0,  are  expanded  in  snapshots,  i.e.  the  columns 
of  matrix  X.  This  can  be  formulated  in  matrix  form  as 


05  =  Xdi,  /  =  !,...,  2m, 


(B  13) 


with  the  column  vector  a,  containing  the  expansion  coefficients. 

The  above  expansion  is  inserted  into  the  large  eigenvalue  problem  (B  12),  which 
results  in  the  2m  x  2m  eigenvalue  problem 

XTWXai=JLciai,  i  =  1,  ...,2m.  (B 14) 


The  eigenvalues  are  the  same  as  the  original  eigenvalue  problem,  and  the  con¬ 
trollable  modes  are  recovered  from  (B  13).  The  orthonormal  set  of  controllable  modes 
are  given  by 


0,  =  —r =Xdi, 

V M 


(B  15) 


There  are  some  important  computational  issues  which  should  now  be  commented 
at:  (i)  The  Gramian  3.2  is  defined  as  an  infinite  integral,  which  means  that  in  order 
for  the  ‘approximate’  Gramian  XXTW  to  be  a  sufficiently  good  approximation,  we 
should  take  snapshots  for  a  long  time.  There  are  no  restrictions  on  how  to  distribute 
the  snapshots  in  time,  and  it  is  prudent  to  store  many  snapshots  when  the  flow  energy 
is  large,  (ii)  Due  to  numerical  round-off  errors,  often  not  all  modes  are  orthogonal. 
In  our  case  with  2/71  =  3200,  the  first  150  modes  were  orthogonal  down  to  numerical 
accuracy  (i.e  (05)r05«  10  15),  whereas  for  higher  modes  the  orthogonality  condition 
was  gradually  lost  due  to  rounding  errors.  The  ratio  /r,  =  k\ //,,■  can  be  used  as  a 
condition  number  of  the  corresponding  mode  0)  .  Very  large  values  of  /x,  indicate 
poor  orthogonality  due  to  numerical  inaccuracy. 

The  observable  modes  are  computed  in  a  similar  manner,  but  now  the  snapshots 
are  taken  from  impulse  responses  of  the  adjoint  equations  for  each  output,  i.e. 
P  =  (Pi,  Pi)  =  )  with  ^  as  the  integrand  in  (2.17).  The  approximate 

observability  Gramian  is 


/»00 

H°i  =  /  fWFVi  dr 
Jo 

where  Y  is  the  n  x  2m  matrix: 

/  px{xx,t i)^  ... 


Y  = 


pT  <b°  dx  dy 


dr 


YYTwj)i, 


Pi(x  1, 


P2(X  1. 


(B  16) 


(B  17) 


\jPl(*^«/2>  •  •  •  P  n/2->  h)^  •  •  •  P  2^n/2->  )  \/^tm  J 

The  observable  modes  are  computed  in  an  analogous  manner  as  the  controllable 
modes  with  Y,  replacing  X  in  (B  1 3)— (B  15). 

B.2.  Snapshot-based  balanced  truncation 

To  obtain  the  balanced  modes,  we  must  diagonalize  the  matrix  SPQ,  which  can  be 
approximated  using  the  matrices  X  and  T,  i.e. 


*  XXTWYYTWtfc  =  CTf0r. 


(B  18) 
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We  expand  the  balanced  modes  as  linear  combinations  of  the  columns  of  X,  with 
a,  =  (ai, . . . ,  am)T  as  the  expansion  coefficients.  Inserting  this  expansion  in  (B  18),  we 
get 

0  =  XXTWYYTWXai-  Xato?  =  X(XT  WYYT  W XUi  -  a, or?) .  (B  19) 

To  solve  this  problem  we  can  equivalently  diagonalize  XrWYYrW X  or  form  the 
singular  value  decomposition  (SVD)  of  YTWX.  The  latter  decomposition  is  preferred, 
since  it  is  numerically  more  stable,  i.e. 

YTWXbi=aiah  1  =  1,...,  2m.  (B20) 

The  normalized  balanced  modes  and  the  associated  adjoint  balanced  modes  are 
recovered  from 

4>?  =  -^=Xbh  tr  =  A=Yalt  (B  21) 

where  (fl°c)T cj)°c  =  Sjh 

The  method  can  be  summarized  in  three  steps:  (i)  Compute  the  impulse  response 
of  each  input;  collect  snapshots  of  the  response;  and  construct  X  (B  11).  (ii)  Compute 
the  adjoint  impulse  response  of  each  output;  collect  snapshots  of  the  response;  and 
construct  Y  (B  17).  (iii)  Form  the  matrix  YTWX;  compute  its  SYD;  and  recover  the 
balanced  modes  from  (B21).  See  Rowley  (2005)  for  further  details  on  the  method. 


Appendix  C.  Timestepper 

The  timestepper  used  in  this  work  for  both  the  forward  and  the  adjoint  equations 
is  a  spectral  code  described  in  detail  in  Chevalier  et  al.  (2007 b). 

If  /( x,  y,  t)  is  a  velocity  component,  then  the  discrete  approximation  is  the 
Chebyshev-Fourier  series 

ny  n*/2 

f(t)  =  Y,T,(y)  Y,  +  c.c.,  (Cl) 

1=0  m=—nx/2 

where  7}  is  the  / th  Chebyshev  polynomial,  am  =  2mn/Lx,  and  nx  =  768  and  ny  =  101 
the  numbers  of  collocation  points  in  each  direction.  The  coefficients  ujm  are  complex 
functions.  The  associated  collocation  grid  is  defined  by  yt  =  (Lv/2)(1 —cos(nl/ny))  and 
xm  =  Lx(l/2  +  m/nx)  with  Lx  =  1000  and  Ly  =  30.  The  discretized  system  of  equations 
is  projected  onto  a  divergence-free  space  by  transforming  to  v-rj  formulation  and 
integrated  in  time  using  a  third-order  semi-implicit  scheme. 

To  retain  periodic  boundary  conditions,  which  is  necessary  for  the  Fourier 
discretization,  a  fringe  region  is  added  at  the  end  of  the  physical  domain  in  which  a 
forcing  is  applied  so  that  the  flow  smoothly  changes  from  the  outflow  velocity  of  the 
physical  domain  to  the  desired  inflow  velocity.  For  the  linearized  equation  the  desired 
inflow  condition  is  zero;  so  the  fringe  forcing  is  of  the  form  F  =  X(x)u,  where 


X(x)  =  —X 


max 


s 


%  %start 
A  rise 


(C  2) 


Here  Amax  is  the  maximum  strength  of  the  damping,  xstart  =  800  to  xenii  =  1000  the 
spatial  extent  of  the  region  in  which  the  damping  function  is  non-zero  and  Arise  =  120 
and  A fan  =  60  the  rise  and  fall  distances  of  the  damping  function.  The  smooth  ‘step’ 
function  S(x)  rises  from  zero  for  negative  x  to  one  for  x  ^  1.  We  have  used  the 
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following  form  for  S,  which  has  the  advantage  of  having  continuous  derivatives  of 
all  orders: 

(  0 ,  x  ^  0 , 

S(x)  =  <J  1/  (1  +  ed/h-d+i/*))  ,  0  <  x  <  1 ,  (C  3) 

1 ,  x>l. 


Appendix  D.  Riccati  equations 

We  briefly  outline  the  full-information  control  and  estimation  problems  and  their 
solutions.  The  reader  is  directed  to  Anderson  &  Moore  (1990),  Lewis  &  Syrmos 
(1995)  and  Bagheri  et  al.  (2008)  for  derivations  of  the  solutions. 

The  first  step  in  the  design  of  an  Jf2  compensator  involves  the  solution  of  an 
optimal  control  state-feedback  problem.  The  full-information  problem  is  to  find  a 
control  u(t)  as  a  linear  function  of  the  flow  state  q(t )  that  minimizes  the  deterministic 
cost  functional 

J  =  \[  qTCTlClq+l1UTU&t  (Dl) 

2  Jo 

while  satisfying  the  initial  value  problem 

q  =  Aq  +  B2u ,  q(t  =  0)  =  q0.  (D  2) 

The  optimal  control  signal  is  given  by 

u(t)  =  -^BlXq(t),  (D  3) 

V - V - ' 

K 

where  X  is  a  solution  of  the  Riccati  equation 

0  =  atX  +  XA-  jiXB2BlX  +  CTlCi.  (D 4) 

The  solution  to  this  equation  provides  the  optimal  steady  feedback  gain  via  the 
relation  (D  3). 

The  second  step  in  the  design  of  an  Jf2  compensator  involves  the  minimization  of 
the  estimation  error  qe  =  q  —  q  given  by  the  estimator 

qe  =  Aqe  +  Blw  +  L(v  —  \z),  q(t  =  0)  =  0,  (D  5) 

v  =  C\q ,  (D  6) 

v  =  Clq+g,  (D7) 


where  w  and  g  are  temporal  white  noise  signals.  The  solution  is  the  feedback  gain  L 
that  minimizes  the  objective  functional 

/*00 

J=  qj(t)qe(t)dt.  (D  8) 

JO 

The  functional  (D  8)  can  be  minimized  if  L  is  chosen  as 

L  =  --2PCt2,  (D  9) 

or 

where  P  is  a  solution  of  the  Riccati  equation 


0  =  AP  +  PAt 


pc2  c2p  +  B\Bl . 


(D10) 
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